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ABSTRACT 

Although the surface of a magnetar is a source of bright thermal X-rays, its spectrum contains substantial 
non-thermal components. The X-ray emission is pulsed, with pulsed fractions that can be as high as ~ 70%. 
Several properties of magnetars indicate the presence of persistent, static currents flowing across the stellar 
surface and closing within the magnetosphere. The charges supporting these currents supply a significant 
optical depth to resonant cyclotron scattering in the 1 - 1 00 ke V band. Here we describe a Monte Carlo approach 
to calculating the redistribution of thermal seed photons in frequency and angle by multiple resonant scattering 
in the magnetosphere. The calculation includes the full angular dependence of the cyclotron scattering cross 
section, the relativistic Doppler effect due to the motion of the charges, and allows for an arbitrary particle 
velocity distribution and magnetic field geometry. We construct synthetic spectra and pulse profiles for arbitrary 
orientations of the spin axis, magnetic axis, and line of sight, using a self-similar, twisted dipole field geometry, 
and assuming that the seed photons are supplied by single-temperature black body emission from the stellar 
surface. Pulse profiles and 1-10 keV spectra typical of AXPs are easily produced by this model, with pulsed 
fractions of ^ 50%. However, this model cannot reproduce the hard, rising energy spectra that are observed 
from SGRs during periods of activity, without overproducing the thermal emission peak. This suggests that the 
1-100 keV emission of SGRs has a common origin with the hard X-ray emission detected from some AXPs 
above ~ 20 keV. 

Subject headings: stars: neutron — stars: magnetic fields — X-rays: stars — radiative transfer — scattering 
— plasmas 



1. INTRODUCTION 

Magnetars have many remarkable properties that appear to derive from the decay of an ultrastrong magnetic field. Measure- 
ments of spindown torques and bright X-ray outbursts have revealed dipo le fields as strong as^ 10'^ G, and imply the existenc e 
of internal fields at least an order of magnitude stronger in some objects (I Thompson & DuncanI 19961 1 Woods & Thompsonl2006h . 
More recent measurements strongly suggest that the external magnetic field can support electric currents that are several orders of 
magnitude larger than the Goldreich- Julian current that flows along the open field lines. In some magnetars, the injection of these 
closed-field currents into the magnetosphere appears to be associated with Soft Gamma Repeater burst activity; but in others, the 
current appears to be sustained in the absence of bright SGR outbursts, in the class of magnetars known as the Anomalous X-ray 
Pulsars. 

The basic properties of a twisted neutron star magnetosphere were outlined by iThompson et alj (l2002h. A detailed one - 
dimensional plasma model of the current in the closed magnetosphere has been constructed bv'Beloborodov & Thompson! (|2006|) . 
The presence of such a global twist has several observational effects, which mirror the behavior of magnetars in the foUowing 
ways. 

1 . Correlation between spindown torque and X-ray hardness. The presence of a twist in the external magnetic field causes 
it to flare out slightly. As a consequence, the open-field flux is larger than estimated from the dipole model, and the spin- 
down torque is increased. The charges flowing along the twisted closed field lines also supply a significant optical depth to 
resonant cyclotron scattering. This optical depth is independent of radius (resonant frequency) when the magnetosphere is 
self-similar. The thermal photons emitted from the surface of the star are multiply scattered. Thus the o bserved correlation 
betwe en spin-down torque and hardness of the power-law component of the 1-10 keV X-ray spectrum (iMarsden & Whit j 
1200 Ih can be explained by adjusting a single parameter: the angle through which the magnetospheric field lines are twisted 
dThompson etal]l2002h . 

2. Persistent changes in X-ray pulse profiles. Following the 1998 giant flare, the pulse profile of SGR 1900+14 changed to 
a single, n early sinusoidal peak, from a more complicated profile with multiple peaks that was observed before the flare 
(IWoods et al. 2001 ). The change in the pulse profile persisted even though the quiescent X-ray flux returned to its pre-burst 
level. Something similar happ ened to SGR 1 806- 20 after its 2004 giant flare, except that instead of simplifying, the pulse 
profile became more complex dWoods et al.l2006l) . The AXP IE 2259+586 also showed a systematic change in pulse profile 
following the 2002 outburst, although the duration of the change was shorter lived than in the SGRs (Woods et al. 2004^ 
It app ears that the external magnetic field of a magnetar is reconfigured following a period of burst activity (IWoods et all 
l200ll:rThompson et al...2002.) . 
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3. Persistent changes in P. The spindown torque of SGR 1900+14 and SG R 1806-20 can show large-amplitude modulations, 
increasing over the previous long term trend by at least a factor ~ 4 (IWoods et alJl2006l and references therein). The 
increase sets in gradually several months after a period of bu rst activity, an d then persists for years. Related but less 
dramatic effects are seen in the torque behavior of some AXPs dGavriil & Ka spi 2004). 

4. Hard X-ray emission. Emission above 20 keV has been measured from both SGRs and AXPs dGotz et aTl l2006l and 
references therein). In the AXPs, the spectrum sometimes shows a strong upward break between 10 and 50 keV, with the 
output at 100 keV being 10 times stronger than the thermal emission from the neutron star surface. It appears that the 
bolometric output of these objects is dominated by the magnetosphere. The additional high-energy component could be 
explained by bremsstrahlung emiss ion from the transition layer near the neutron star surface, or possibly by a pair cascade 
at ~ 100 km from the neutron star ([Thompson & Belobor odov 2005). 

In this paper, we confront the magnetar model with some of the observational evidence listed above. We focus on the part of 
the non-thermal X-ray spectrum that lies directly above the black body peak, and also on the X-ray pulse behavior in this band 
(items 1 and 2 above). The most robust way of testing the twisted mag netosphere hypoth esis is to produce synthetic analogs 
of observable quantities, using the self-similar model developed by Tho mpson et al.l (l2002h . The radiation transport through the 
magnetosphere is sensitive to several geometrical effects, which include the angular emission pattern of cyclotron scattering by 
moving charges, the magnetic field geometry, orientation with respect to the observer, and the particle velocity distribution. 

A first strategy is to ignore all the geometrical complications, focusing on the essential physics entering the multiple scattering 
process in a one-dimensional geometry. The problem of radiation transport through the accretion colu mn of an X -ray pulsar has 
been treated in this way, using both Feautrier methods to solve the r adiation transport equation (Mesza ros & Nag el 1985a b) and 
Monte Carlo techniques ( Pravdo & Bussardlll981l : IWan g et al. 1988). The more general problem of the transfer of X-ray photons 
through a single resonant surface has been analyzed by Zheleznyakov ( 1996) and Lvutikov & Gavriil (2006) by restricting the 
transfer of photons to one spatial dimension and assuming a sub-relativistic and thermal distribution of scattering charges. 

We do not expect narrow spectral features to form as the result of the transfer of thermal X-ray photons through the quiescent, 
twisted magnetosphere. The cyclotron optical depth is regulated to a value close to unity over a broad range of frequencies. By 
contrast, the optical depth is much higher in a pulsar accretion column, and is concentrated over a narrower range of radius and 
magnetic field. (In addition, the equivalent width of the ab sorption fe atur e that forms in the cold atmosphere of a passively cooling 
magnetar is suppressed by polarization mode exchange iHo & Laill2004i ) The optical depth varies significantly with angle in the 
twisted dipole magnetic field, and the field lines impart a non-trivial angular structure to the particle velocity distribution. This 
means that a numerical simulation in three dimensions is required to model the formation of a high-energy tail to the spectrum, 
and the amplitude and shape of the X-ray pulses that are seen by a distant observer. 

In this paper, we report the development of a Monte Carlo code that follows the transfer of X-ray photons through the stellar 
magnetosphere, with resonant cyclotron scattering as the source of opacity. The code treats the scattering problem in three 
dimensions, for an arbitrary magnetic field geometry, and an arbitrary velocity distribution of the scattering charges. We account 
for the relativistic Doppler effect and polarization exchange during scattering, but neglect the effects of recoil (which are important 
only at high energies where the spectra of AXPs and SGRs are not well measured, and where additional emission components 
are sometimes observed). The output spectra and pulse profiles are obtained for an arbitrary line of sight, and for an arbitrary 
orientation of the spin axis and magnetic axis of the neutron star. 

The paper is structured as follows. In § 2 we review the physical input to the model, and in § 3 we describe the Monte Carlo 
code. Results are presented and discussed in detail in § 4. A comparison with the existing data and general conclusions follow in 
§5. 

2. RESONANT CYCLOTRON SCATTERING IN A NON-POTENTIAL MAGNETOSPHERE 

In this section we describe the physical model that w e are simulating. We first review the properties of a neutron star magne- 
tosphere that carries a net twist ([Thompson et al.ll2002l) . We then review the process of resonant cyclotron scattering and how 
it depends on the structure of the magnetosphere. At frequencies lower than the surface cyclotron frequency, the optical depth 
depends only on the twist angle and the drift speed of the charges. The effects of relativistic particle motion will be considered in 
detail. 

In the magnetar model, there are two sources of keV photons. The first is deep cooling. As the internal magnetic field is 
transported outwards, part of the field energy is converted to heat, and the remainder to a deformatio n of the magnetic field outsid e 
the star Part of this internal heat is conducted to the surface, and the remainder radiated as neutrinos ("Thomp son & Duncanll9 96h. 
This component of the blackbody flux will be carried by the E-mode, given that its opacity is strongly suppressed in the surface 
layers by a factor ^ (meCUj/eB)^ ^ (^B7bb/»^^.c^)^B75 (iSilantev & Iakovlevli980 ). 

The impact of curr ent-carrying particles on the surface is a second possible source of kilo-electron volt photons 
([Thompson et al.ll2002h . Some fraction of the particle kinetic energy will be deposited in a layer which is optically thick to 
the O-mode, but optically thin to the E-mode. This energy would then be re-radiated in the O-mode with a nearly black body 
frequency distribution. Although ions are probably able to penetrate this deep, it has been argued that a relativistic electron (or 
positron) beam will deposit a signficant frac tion of its energy in an optically thin layer, which is re-emitted in the O-mode at 
a much higher temperature ( Thompson & BeToborodovl[2005h . A recent re-cali bration of the distance to sever al AXPs shows 
that their 1-10 keV luminosities are tightly clustered around ~ 1 x 10"^^ ergs s"' ([Durant & van Kerkwiik[|2006bh . These sources 
have a range of spin-down rates and their X-ray spectra are of variable hardness. The clustering of the 1-10 keV emission 
therefore suggests that the rate of heat conduction to the surface is buffered by bulk neutrino cooling, and that the black-body 
emission is not dominated by particle heating. In at least one AXP, the magnetospheric output at ^ 100 keV is significantly 
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bri ghter than the ~ 1 ke V surfac e emission (iKuiper et al.ll2005Ll2006h . and is more comparable to the expected neutrino luminos- 
ity ([Thompson & Duncanll 19961) . This suggests that a significant fraction of the energy that is released by the untwisting of the 
internal magnetic field ends up in a quasi-static deformation of the external magnetic field, such as we are assuming in this paper. 

A photon propagating away from the surface of a magnetar can scatter resonantly at a significant rate over some range of its 
path. If the magnetospheric particles had a narrow velocity distribution of width A/3c, then the cyclotron scattering would be 
localized on a 'resonant surface' with a well defined radius R^es and a thickness A/?ies ~ A/3/?ie,s. A similar statement holds for 
any photon after it has been resonantly scattered within the magnetosphere. 

The transfer of thermal photons through a single warm scattering layer (with one-dimensional electron temperature 
^bT \\ <C mgC^) results in a second component of the spectrum that is upshifted in frequency by a factor 1 + lik^T^^ / m^c^y 
(iLyutikov & Gavriiill2006i) . This analytic result is, however, restricted to one spatial dimension, does not include bulk motion of 
the charges, neglects the effect of scattering at distant points in the magnetosphere, and becomes inaccurate when the r.m.s. speed 
of the charges exceeds^ f3 ~ 0.3. 

In a more realistic, three-dimensional magnetosphere a power-law component of the spectrum results from two effects: 

1. Multiple resonant scattering within a single magnetic hemisphere, where the charge flow is converging and the average 
increase Auj in photon frequency per scattering is first order in (3; and 

2. Multiple scatterings in different magnetic hemispheres, for which Aw ~ f3^ if the particle drift speed parallel to the magnetic 
axis is symmetric in the two hemispheres. 

For example, when the charge flow in the twisted magnetosphere is carried by electrons and positrons, each hemisphere will 
contain one species of charge whose bulk motion is directed toward the star Even when electrons are the only species of light 
charge (the positive charge flow being carried by ions), the electron velocity will converge in one magnetic hemisphere. The 
optical depth is independent of frequency in a self-similar magnetosphere, and a significant fraction of the backscattered photons 
undergo multiple scattering if the field is twisted through ^ 1 radian (Thompson et al. 2002; eq. ifTTIl '). 

2.1. Twisted Self-Similar Magnetospheres 

[Thompson et al J (|2002|) solved the force -free equation J x B = outside a spherical surface with a self-similar ansatz 
(iLvnden-BeU & Boilvlll994 IWolfson|[T99l . These solutions form a one-parameter sequence, labeled by the net twist angle 
of the field lines that are anchored close to the magnetic poles, 

A(/)N-s=2 1im / -^TTTK^-R- (1) 



0^0 Bg{6) sin 6* 

Here 9 and (p denote the polar and azimuthal angle relative to the magnetic axis. The magnetic field is 

B(r,0) = ^(^y''F(cos0), (2) 

where Bpoie is the field strength at the magnetic poles, r is the radial coordinate, and /?ns the stellar radius. The components of F 
are expressed in terms of a function /(cos 6), which is the solution to the ordinary differential equation 

^m^ef" + Cf'-^l" + p(p+\)f=0. (3) 

One finds /v = -/'; Fg = {pj sin6')/; and = [C/p{p+ ^ Fo- The function f(cos9) satisfies the three boundary 

conditions /'(O) = 0, /'(I) = -2, and /(I) = 0. These force-free equilibria^ smoothly interpolate between a dipole (A0n-s = 0, 
p= 1) and a split monopole (A0n-s = tt, /? = 0). The poloidal field lines remain close to a dipole as long as A0n-s ^ 1- 
The current density induced along the twisted field lines is ([Thompson et alj|2002h 

J=\Zienif3iC=— -^B~— sm 6IB, (4) 

^-^ 4-TTr Bg Airr 

i 

where Z,e is the electric charge of particle species /, n, is the particle density, and (3^ = (3iB is the particle velocity directed along 
B = B/Z? in units of the speed of light c. 

The closed-field current ^ implies a much higher flux of charged particles than can be supplied by the rotationally- 
induced (corotation) ch arge density. As a result, the ma gnetospheric plasma contains both positive and negative charges with 
nearly equal densities (iBeloborodov & Thompson! l2006l) . One has y/|pGj|c ^ A0n-s (f^^max/c)"', where pcj = •B/27rc 
dGoldreich & Julianll 19691) . fl is the spin angular velocity of the star, and R^ix — r/ sin^ 6 is the maximum radius of the poloidal 
field line that passes through the position (r, 9) in the magnetosphere. We will consider both a unidirectional flow of the scattering 
charges at any given position in the magnetosphere, and also a bi-directional flow. The first case corresponds to an electron-ion 
plasma, and the second to a pair-dominated plasma. 

Working to first order in /3, one must approximate [1 + (<:B7'||/meC-)'/^]/[l — (AibT"!! /"icC^)'^^] — 1 +2(i:Br|| /nipC^)'''^. 
** Tliere are two solutions for each value of the constant C < .873, and the two branches connect smoothly at C = .873. Together they form a one-parameter 
sequence labeled by p or A<jiN_s- 
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2.2. Basic Properties of Resonant Cyclotron Scattering of X-ray Photons 
When a plasma is threaded by a magnetic field, cyclotron absorption will occur whenever the condition 

LU = UJc— (5) 

mc 

is satisfied in the rest frame of a charged particle. Here, oj is the photon frequency and Ze and m are the particle charge and mass. 
Absorption by motionless charges occurs at a radius 

1 l/(2+p) 

(6) 



/3 
■ole, 

the angular dependence of the magnetic field strength: 



1/3 1/3 

For a 1 keV photon, this works out to /?res ^ ^^^^poie 15 electrons and 205p^jg ^3 km for protons. The function ^ contains 



e = -^fi^ (7) 



l + 3cos2 6' 



for a dipole, and 



^ 4 



smO J p+1 



(8) 



for the self-similar twisted magnetosphere (eq. IS))- 

The lifetime of the first Landau level, 1/r = 3OT^cV(4e'*B^) 3 x 10"'*(OT/m(,)^(B/10'^G) s, is much shorter than the dy- 
namical time Rres/c when the cyclotron energy is greater than ^ 10~^(m/me)'*''^ eV. In this fast-cooling regime, the combined 
absorption and reemission can be treated as a scattering process, with a rest-frame cross section 



(e.g. lMeszaroslll992h . The amplitude 



fTres=4^'^k/..pt^.<5(w-Wc) (9) 



e,,, = e-mf] (10) 



V2 



is the overlap of the photon's polarization with a left (+) or right (-) circularly polarized mode (depending on whether the 
scattering charge is negative or positive). Here i and y are unit vectors perpendicular to B. Integrating through the cyclotron 
resonance yields a cross section that is much larger than Thomson, by a factor ^ {e^ /Tic)~^{huOc/meC^)~'^ ~ 10^ for a keV photon. 
We neglect the electron recoil in our calculations, and relativistic effects on the scattering cross section (e.g. iHeroldl 19791) . These 
effects start to be significant at photon energies > 50 keV ^ mgC^/lO if the electrons are sub-relativistic. 

The surface magnetic field of the AXPs and SGRs is believed to exceed Bqed = m^c^ /eh = 4.4 x 10'^ G, where the energy 
of the first electron Landau level begins to exceed ^ nieC^. We are concerned here only with resonant scattering in the parts of 
the magnetosphere where the Landau levels of the scattering charges are non-relativistic. Resonant absorption of kilo-electron 
volt photons close to the surface of a magneta r creates a high-energy gamma ray that is converted to an e"*" - e~ pair near the 
emission site dBeloborodov & Thompsonl2006h . The ion cyclotron energy is in the kilo-electron volt range even at the surface of 
a magnetar, and so the effects of ion cyclotron scattering are concentrated much closer to the star The cross sections that follow 
generally apply only in this non-relativistic regime. 

The current ^ supplies enough particles for resonant cyclotron scattering to become the dominant source of opacity in the 
keV range (.Thompson et al..,200Z; .Lvutikov & Gav riil 2006). Integrating eq. (|9]l through the resonance gives an optical depth 

A0N-S . 2q .ns 
Tres 1^ Sm 9. (11) 

This expression does not depend on the mass or charge of the scattering particle, and applies as long as the resonant energy is less 
than mc^ and lies below the surface cyclotron energy. This result is a dkect consequence of the self-similarity of the force-free 
construction. A more realistic magnetosphere has a more inhomogeneous distribution of twist, contains higher-order multipoles, 
and is non-axisymmetric. 

The resonance condition is changed by bulk motion of the charges along the magnetic field fines, and becomes 

^ = ^^0^—. — (12) 

Here 7 = ( 1 - /3^)"' is the Lorentz factor and ii = k-B = cos 0kB is the direction cosine of the photon with respect to the magnetic 
field in the stellar frame. In this case the resonant radius in eq. (|6]l needs to be m ultiplied by a factor [7(1 -/ ?/i)1~^/^^"'"^'\ Similarly, 
the cross section in the stellar frame equals eq. (|9]l multiplied by ( I — f3^) [e.g. iRybicki & Lightrnanlll979ll . with ujc replaced by 
wd, and with |e;., P evaluated in the rest frame of the charge: 

a,,, =471^1- (3 ^^)^-^\e',^\^LUuS(Lo-LOu) (/?^0). (13) 



Quiescent X-ray Emission of Magnetars 



5 



Here the ' labels the particle rest frame. 

The polarization modes of a photon i n a magnetized pl asma depend on the relative contributions of plasma and vacuum 
polarization to the dielectric tensor (e.g. IZheleznvakovlll99 6). Even in the core of a cyclotron resonance, vacuum polarization 
dominates in a cold plasma with a velocity width A/3 <C 1 if the condition 



LOT: 



(B < Bqed) 



(14) 



is satisfied. In eq. (fT4b . ujp = ^ A-Ke^nejme is the electron plasma frequency. This means that both electromagnetic eigenmodes 
are then linearly polarized even at frequencies very close to lod. The extraordinary mode (E-mode) corresponds to e = B x ^, and 
the ordinary mode (O-mode) to e = (B x k) x k. This means that each polarization eigenmode will experience the same interaction 
with the gyromotion of a positron as it does with the gyromotion of an electron: the overlap 



\e, Al- 



l/2; =cos2 0kB/2 



(15) 



does not depend on the sense of circular polarization. In the case of unpolarized photons, 



.|.,„„„l^(l+ C0S^gl,B)/4. 

Mode exchange is clearly possible during resonant scattering (e.g. IZheleznvakovll 1 996t IWang et al]|1988l) . The differential 
cross section for scattering from mode A to mode B is dap^^/cKl' oc |ej ^(A)p|eJ ^.(B)p. 
converted to the extraordinary mode is 



The probability that any mode A is 



P{A E) 



^;..(E)p 



1 



and the probability for conversion to the O-mode is P{A 
scattering) between k and B in the charge rest frame. 



<,(E)|2+|<,(0)|2 

0)= l-P(A- 



l+cos2 61 



(16) 



E). These branching ratios depend on the angle (after 



(17) 



2.3. Basic Degrees of Freedom of the Magnetospheric Model 



The model is defined by 



1. The twist angle A(/)n-s (eqs. ID, ifTTl ). which labels a one-parameter sequence of magnetic field geometries. We con- 
sider only axisymmetric field configurations in this paper (Our numerical implementation can in principle handle more 
complicated, non-axisymmetric, geometries.) 

2. The spectral distribution of the seed photons. We generally assume a black body distribution, defined by a single tempera- 
ture 7bb. In numerical calculations, smoother spectra below the black body peak are obtained by inputing photons at a single 
frequency, and then convolving the output spectrum with an input black body spectrum. We show in Appendix |A] that this 
procedure nicely reproduces the spectrum that is obtained from drawing photons dkectly from a black body distribution. 

3. The polarization of the seed photons. We calculate the output spectrum assuming that the seed photons are 100% polarized 
(E-mode or O-mode). The spectrum resulting from an arbitrary admixture of seed polarizations can be obtained by linear 
superposition. 

4. The angular distribution of the seed photons. We assume that the emitting surface is homogeneous. The output pulse pro- 
files and pulsed fractions depend more sensitively on the surface temperature distribution than does the high-energy tail of 
the X-ray spectrum. We wish to define the pulse profile shapes that can result from radiative transport in a current-carrying 
magnetosphere, and the maximum pulsed fractions that can be obtained without resort to temperature inhomogeneities 
across the surface. 



The ratio of the seed photon energy to the cyclotron energy of the scattering charges at the magnetic pole, e.g. 
^B7bb(^|Ze|Z?poie/'«c)"'. Although |Ze|, Bpoie, and m are all independent quantities, they enter only in this ratio. One 
must, equivalently, specify the ratio of /?ns to rbb = ^res(w = ^B^bb/^) (eq. IS)), 



^ = (7x10-2)'/^'"'^' 



1 



feB^bb \ / lO'^G 
0.4keVy Bpoie 



l/(2+p) 



(18) 



Because this ratio is small for electrons and positrons, in our calculation we inject photons moving in the radial direction 
from the surface of the star. The propagation of a non-radially moving photon would result in nearly radial propagation 
at the point of first scattering. We have checked (directly by simulation) that injecting the photons with a broad angular 
distributions results in practically identic al sp ectra between ^ 1 and ^ 50 keV. A broader range of injection angles is also 
considered in the case of ion scattering ( ^4. 5b . where it is shown that there are modest differences in the output spectrum. 

6. The momentum distribution f{p) of the charge carriers. To explore the dependence of the results on /(/?), we use three 
distribution functions, which are taken to be independent of position: 
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I A mildly relativistic flow in which the particle motion is restricted to one direction along the magnetic field (e.g. 
j3 > 0). We choose a one-dimensional Boltzmann distribution in energy, which is labelled by a temperature kBTo = 
(70- 1)otj,c^, or equivalently by a drift speed (3o = (1-7(7^)'/^, 



= M V n /r TT: 

A^/3^i(l/[7o-l]) 



7o-l 



A^;3^i(l/[7()-l]) 



exp 



7o-l 



(19) 



Here Ki is the modified Bessel function, and Np number of directions in which the scattering particles move along 
the magnetic field. This is a distribution in momentum /37, and is normalized in the usual manner; 

d{[3-i)f{(}-i)=\ (20) 

Jo 

for the asymmetric distribution with Np = \. 

II A mildly relativistic, one-dimensional gas with the same Boltzmann distribution as I, but now extending over positive 
and negative momenta, -00 < (3^ <oo. In this case, the normalization of / is smaller by a factor of A^^' = j . This type 
of symmetric distribution plausibly describes an electron-positron gas. In the magnetospheric circuit, the positrons 
are identified with the half of the distribution function with positive velocity, and the electrons with the particles with 
negative velocity. (Both types of particle must exist in essentially equal numbers when the plasma is dominated by 
pairs, because it is nearly charge-neutral.) 

Ill A broad power-law in momentum, 

/(/37) « (/37)", (21) 

with minimum and maximum values {fi^)mm and (/37)max, representing a warm relativistic plasma. The index a = -1 
corresponds to an equal number of particles per decade of momentum. We consider only distributions with a < 0, 
since distributions with positive a are subject to vigorous kinetic instabilities (e.g. Kulsrud 2005). 

These three families of distribution functions are characterized by the parameters {/?()} for the Boltzmann distributions, 
and by {/3min,7max, ct] for the power-law distribution. In practice, the distribution functions are sampled only over a finite 
range of velocities, and upper and lower cutoffs are applied to the admissible values of /3. For the distributions of type I 
and II, these cutoffs are taken to be very close to zero and to unity, respectively. See 33.4l for further details. 

The particle distribution function must actually vary with position in the magn etosphere. A strong drag force act s on electrons 
and positrons at r ^ 100 km, where their cyclotron energy is in the keV range ( [Thompson & Beloborodov 2005*). As a result, 
one expects the motion of the light charges to be mildly relat ivistic in this part of the magneto sphere. Published calculations 
of the charge dynamics in the closed magnetospheric circuit dBeloborodov & Thompsonll2006h allow for pair creation but do 
not model the outer part where B < Bqed- Close to the star, the particles have a broad relativistic distribution, extending up to 
an energy 7ies'WeC^ ^ (TjeBpoie/'WeC^BTbb) dBeloborodov & Thompsonll2006h . At this Lorentz factor, it is possible to spawn new 
electron-positron pairs through the excitation of a charge to its first Landau state, by absorbing a photon from the peak of the 
thermal radiation field. Understanding the interplay between these two competing effects require coupled radiation transfer and 
electrodynamic calculations, which we will pursue elsewhere. 

3. MONTE CARLO SIMULATION 

The random nature of photon scattering in the magnetosphere, and the statisti cal character of observable quantities (spectra, 
pulse profiles) makes the Monte Carlo technique especially suitable (see, e.g., ISoborllT994l for an overview of Monte Carlo 
methods for Compton scattering in hot plasmas). We now present the details of the magnetospheric model and the numerical 
algorithm. The results of our calculations are reported in ^ Some validation tests are described in Appendix lAl 

3.1. Transfer of Photons through a Cyclotron Scattering Atmosphere 
We now present the basic equations that describe the transfer of radiation through particles with an arbitrary (relativistic) 
distribution function. A photon following a ray path r = ro + £k sees a differential optical depth 

(ifi-es(w,fe,r) f\ dnz - 

d(3 — a,^,(uj,k,r) (22) 



di J_i d(3 

for cyclotron scattering. The point of emission of the photon is labelled by ro (t = 0). The one-dimensional velocity j3B of the 
charges has the distribution 

/z(Ar) = «z'^, (23) 

where j'j fz(f3,r)dP = 1. Given that particles of charge Ze carry a fraction ez of the current, their space density nz can be 
obtained from eq. (|4|i 

\Z]enz (P+IXB, 
B UnrlPl Be ' ^ ' 
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We focus on a single species of scattering charge in the remainder of this section, and drop the subscript Z. Radiation transfer 
through an electron-positron plasma is easily described by writing n = n-+n+ and f{(3,r) = /-(/3,r)+/+(/3,r), where the subscript 
labels the sign of the charge. We will treat the effects of ions and electrons separately in the case of an electron-ion plasma. In 
the formulae below, the factor e represents the fraction of the current that is carried by the scattering charges: e = 1 for a pure 
pair plasma. Substituting expressions (fTSl l. (|23] |. and (l24l i into eq. (l22l l gives 



dr. 



_res ^ 7r(/:»+l) fB_ 
dl ~^ r \B, 




(i-/3M)k,'J 



' |2 



■uoo6{uj-ujv))fifi,r)dp\ 



At a given location in the magnetosphere, the resonance condition (fT2] i can be satisfied only for two values of (3 



(25) 



(26) 



which depend on the photon frequency uj and propagation direction k, and the local cyclotron frequency uJc- We can eliminate 
the delta function in dZSl l through the substitution 



6{uj-ujb) = 



because uju is function of f3 only for fixed r, uj, and k. The velocity-averaged optical depth is then 



dTi-sfi Tr(p+l)uj 



d£ 



= e- 



(1-/3V)|<J2^ /(/3^r) (l-/3-M)|<J|- /(/3-,r) 



\diOu/d(3\p. 



1/3- 



(27) 



(28) 



The superscripts correspond to the two roots (3^ . 

Once the frequency and direction of the photon are fixed, the resonance condition can be expressed in terms of two functions 
(3^ of two local parameters lOc/uj and fi. A narrow distribution of charges in velocity space implies the existence of a well-defined 
'resonant surface' in coordinate space (Fig.[T]i. (The shape of this resonant surface shifts with the direction of the photon if the 
charges are in bulk motion.) Resonant scatterings are possible only in a restricted part of this parameter space: eq. (|26] | has 
solutions only if 

M^+(^J)'>1. (29) 

The division of the photon trajectory into discrete steps, whose size is adapted according to a pre-determined velocity distribu- 
tion function, is the core of the Monte Carlo method used in this paper The propagation of a photon through the magnetosphere is 
calculated most conveniently in this resonant velocity space: this allows the step size to be weighted appropriately by the density 
of the resonating charges. See the left panel of Fig.[T]for an illustration. We neglect the effect of gravity on the photon trajectory, 
so that k = constant in between scatterings. Applying the resonance condition ijJo = 'jJ (eq. ifTSlD gives 



duJr 



cjd dB dujD dfi 
'Y'dl^^'dl 



0, 



(30) 



where dBjdt = V5 • k and dfi/dt = V/x • k. The relation between the step size in coordinate space and in velocity space becomes 

(5c.d/5/3)^± 



M = - 



where 

The differential optical depth is then 

Tr(p+ l)u! ( 



1 duj-D 1 dB 



UJ-D 



d£ 



B 



13 d^ 



ATres = £- 



B, 



(i-rM)k;..ii /(r,r) 



(31) 



(32) 



(33) 



m {duJo/dl)^, " 1/3-1 {dLO^/di)p_ 

Since duj^ / d£ is proportional to uju = the optical depth is independent of frequency. The two terms in (l33T l have opposite signs 
because the two branches of (3^ advance in opposite directions in velocity space for a given photon direction in coordinate space 
(see left panel of Fig. [Hi. When substituting eq. dSTl i into eq. ( |28] |. we get a factor 

V|9c.d/9/3|L± 



(34) 

) is generally dominated by 



which is +1 for j3~ and -1 for (3'^. The gradient of the Doppler-shifted cyclotron frequency (eq 
the gradient in magnetic field strength, and is negative for photons moving away from the star. 

The photon step size A£ can be chosen in such a way that the products r)A/3^ (as determined from eq. 13TI ) yield a 

small value for Afies . This method ensures that the bulk of the resonant surface (where /(/3^ , r) is high) is traversed with a small 
step size. 
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Fig. 1 . — Left panel. Solid curve denotes the velocities (eq. 1261 ) of the charges that are in resonance with a photon of frequency and direction cosine 
/i (as measured with respect to local magnetic field). Vertical dashed line represents the escape surface, which is determined self-consistently by the condition 
u)cl'^ = \/ 1 — /.t- , or equivalently by = (see text). No resonant scattering is possible to the left of this line. When the velocity distribution of the charges 
is narrow (horizontal shaded bars), resonant scattering is possible only with a narrow range of i^d^ and of radius. This resonant surface is even narrower in 
coordinate and frequency space when it overlaps with the escape surface. Right panel. Sample photon trajectory in the {ll)i:/u))—ii plane. Here the escape surface 
is represented by the dashed curve + (uJc/lo)^ = 1. The photon trajectory starts at the surface of the star (to the extreme right, not shown), enters the plotted 
domain at point A, and scatters at point B. After the first scattering, the frequency and direction cosine of the photon shift over to point C. In the example shown, 
the photon is backscattered to smaller radius, and leaves and re-enters the plotted domain at points D and E. It finally escapes the magnetosphere at point F 
without any further scatterings. 

3.2. Photon Escape from the Magnetosphere 

A photon cannot resonantly scatter in the outer region where p^ + {u)c/lo)^ < 1 (see eqs. Il26l and ll29l ). This inequahty is 
a sufficient condition, which is independent of the velocity distribution of the charges. Nonetheless, it is possible for a photon 
to leave and then re-enter the inner region where resonant scattering is allowed. To see this, it is easiest to examine the photon 
trajectories in the {uJc/lo)- p plane (the right panel of Fig. [U. The trajectories are generally not straight lines in this plane, but 
have a mild curvature close to the escape surface, as long as no scattering occurs. 

We therefore adopt a second criterion, which must be satisfied if a photon in a Monte Carlo simulation can be considered to 
have entirely escaped the star and become observable. The tangent vector of the photon trajectory in the {ud lo)- p plane 



t: 



ujc \ 1 dB . dpi , 

ijJ-\ u 

uj J B dl de^ 



(35) 



must not point in a direction that intersects the escape surface for a second time. [Here the unit vector uj runs parallel to the 
(ujc/u!)-a.xis.] In other words, we require that t point in the direction of the line segment ipjc/uj) = 0, -1 < /i < 1 that forms the 
left boundary of the (ajf/w)-/! plane. If this condition is not satisfied in combination with /i^-l-(cjc/aj)^ < 1, then it is necessary 
to follow the photon trajectory to check whether it re-enters the region where resonant scattering is allowed. 

A special case occurs for monopolar fields, or photons emitted very near the magnetic poles in a dipole field. In both these 
cases, the photon trajectory will tend to /j, = ±1 at large distances from the star. In numerical calculations an arbitrary cutoff 
needs to be taken. In our model this is not a problem for dipole-like fields, since the optical depth vanishes at the magnetic poles 
(recall eq. E)). 

The curvature of the resonant surface plays an important role in determining the rate of photon escape. Consider the case where 
~ p, so that the resonant surface of a photon overlaps its escape surface. For the purpose of illustration, imagine a top-hat 
velocity distribution that is independent of position, and has a narrow width A/3 centered about (3 = (3, 

1. \f (R- AR l'->\ ^ R ^ (Ra- \R 

(36) 



/(/3,/3,A/3) = 



^ if (/3-A/3/2X /?<(/? + A/3/2); 
otherwise. 



Two such distribution functions are depicted as horizontal shaded bars in the left panel of Fig.[T] If (3 happens to be very different 
from /i (the lower example), then the resonant surface has a width 



A/3 



d\nB 



d\nr 



(37) 



which depends on A/?. But if (3- i A/3 < p < /3+ 5A/3, then the resonant surface becomes very thin in the radial direction, and 
relatively thick in the non-radial direction. The distance over which the photon escapes from the resonant layer then depends on 
the rate at which the layer curves away from the photon trajectory. 

In practice, enough photons are scattered into an angle /i — /3 that the transport of photons through a resonant layer is not well 
approximated by a one-dim ensional mod el (e.g. a model in which the cyclotron frequency has a gradient only perpendicular to the 
resonant surface: IZheleznvakovul996;iLvutikov & Gavriilll2006i) . In this situation, the step size has to be reduced appropriately 
(see next section). 
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Fig. 2. — Flow chart showing evolution of single photon. Ellipses contain starting, ending, and loop points in the algorithm, diamonds contain decision points, 
and rectangles contain simple instructions. See text for a detailed explanation. 

3.3. Code Description: Single Photon Evolution 

We evolve a single photon at a time, from the point of emission on the neutron star surface to the final escape (as defined by 
the criteria § I3.2| |. The final frequency Wout and direction cosine /i^ relative to the magnetic axis are then recorded. This process 
is repeated until the distributions in Wout are are well sampled, according to criteria discussed in Appendix A. 

In our code, a photon is a seven-dimensional vector, the components of which are its position r, direction k, frequency lo and 
polarization (E or O).^ Although the photon's electric vector is defined with respect to the local magnetic field, the polarization 
will adjust adiabatically as the photon propagates through the magnetosphere.* The code is implemented in FORTRAN??. 

Each photon is injected at a random point on the stellar surface, r = Rj^s, with a wavevector pointing in the radial direction, with 
a frequency drawn from a black body distribution, and with fixed polarization. Once the initial ray path is chosen, the photon's 
position and frequency are updated through a set of binary decisions: see the flowchart in Fig.|2] Except for extreme values of 
the initial frequency, the scattering rate of the photon is very small just above the point of injection, nevertheless, we assume very 
small seed photon step size (/?ns/50). 

First, the subsequent step size along the photon's trajectory is determined by an adaptive algorithm ( ^3. 4b that minimizes the 
step size in regions where the density of resonating charges is highest. The step size is further reduced whenever ~ fi, until 
the condition 

< 10- (38) 

is satisfied. 

Second, we check whether the escape criteria defined in § 13.21 are satisifed. If the photon must escape, its frequency ui and 
direction angles /ik and (fn^ are recorded. (The magnetic field is axisymmetric in the calculations presented in this paper, and (/)k 
is not used in the analysis.) 

If the photon does not escape immediately, then we determine whether it scatters by drawing a random number (Appendix |3.4t . 
The contributions of the two resonant populations of charges (3^ are included in the calculation of the differential optical depth. 
If the photon scatters, then its vector is updated. A new value for the step size is calculated, the photon is advanced, and the loop 
is restarted. 

If the photon does not scatter, we advance it once step. Before restarting the loop, we check whether it hits the stellar surface. 
An E-mode photon that hits the surface is assumed to be absorbed: its non-resonant opacity is strongly suppressed, and so the 
photon penetrates to a depth where the absorption opacity dominates over scattering. Such a photon does not contribute to the 
output spectrum, and only its absorption is recorded. We allow half of the O-mode photons hitting the star to scatter elastically 
off the surface, and half to be absorbed. 

Each photon is assumed to propagate on a straight line, and the effect of gravitational light bending close to the star is ignored. 
This approximation is excellent near the first resonant scattering surface if the scattering charge is an electron or positron, since 
Tres ^ 2GMns/c^ ^ 27?ns/3 (scc eq. ||6l); here Mns is the neutron star mass). The temperature is assumed to be constant over 
the stellar surface, and so the effects of light bending can be neglected closer to the star Light bending is more important when 
the ions are the dominant source of opacity in the magnetosphere (as they may be in the period of strong afterglow following a 

^ The photon's position and direction are defined in a spherical coordinate system centered on the star, with z parallel to the star's magnetic axis: r = (r, cos 9 = 
f-z, (/>) and k = (/i^ = k-Z, rpk)- Our unit of frequency is aj),b = ^B^bb/^. of magnetic field Bbb = nicuj],i,/{\Z\e), and of length rbb = Rres('^bb) (si- (HI)- 
* The polarization tracking of an X-ray photon breaks down only outside the last surface of cyclotron scattering. 
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Fig. 3. — Angles defined at a scattering site, in the stellar frame. Tlie vectors z, B, and k coirespond to the directions of the magnetic axis, the local magnetic 
field, and the photon wavevector, respectively. See text for the definition of other quantities. 

Soft Gamma Repeater Bursts: Ibrahim et ani2001h . We present some preliminary results for ion scattering in ^4.51 The combined 
effect of light bending and surface temperature variations will be included in future calculations. 

This Monte Carlo procedure is independent on the magnetic field geometry and velocity distribution function, as long as only 
one species of scattering particle is involved. It can easily be generalized to include multiple scattering species, and multi-peaked 
velocity distribution functions (electron-positron pairs). 



3.4. Adaptive Step Size Algorithm 

The ev olutio n of a photon's position is performed in velocity space for type I and 11 distributions, and in momentum space for 
type III 02.31 ). To maximize accuracy and minimize computing time, larger steps are taken where /(/3*,r) is small, and vice 
versa. 

To implement this approach, we assume that any velocity distribution at a given position r can be characterized by a single 
peak with half-widths A/3iow and A/3high, centered at /? = In the case of a broad distribution, we take the median value of [3 as 
the central point. In the case where distribution function is supported over a range of (3 extending downward to or upward to 
unity, (e.g. the Boltzmann distribution), we restrict the values of (3 at the extremes of this range to avoid numerical problems. The 
limits are chosen empirically so that the output spectra and pulse profiles are insensitive to the choice. In practice, this involves 

sampling some 99.8% of the (integrated) distribution function, so that J^"" f(P)dl3 = Jp^^ J(l3)dP = 0.001. Thus, if both /3+ 

and /3" are outside the resonant region, we take a big step in coordinate space (A^ = r/10). If the value of either /J""" or /3" is 
approaching a possible resonance, we refine the step in velocity/momentum space until the edge of the region is resolved to a 
fiducial accuracy, namely APlg^/^^[^^^/lQQ. Finally, if either of them lies within the resonance surface, we set the step to a small 
logarithmic increment in momentum space: 

dhf3\,,^ = sgn[flf(7/3)oid] hP\/lOO, (39) 

or equivalently, df3new = sgn((i/3oid)(l-/3^)|/3|/100, where the subscripts old and new refer to the stepsize before and after refining, 
and sgn is the sign function. The numerical factor was calibrated empirically to ensure reasonable smoothness in the angular 
distributions of all cases explored. For the special case (3^ = 0, we set the step to a small constant. Once the appropriate size of 
the step in either coordinate or velocity space has been determined, the complete set of steps (A£,A(3'^, and Af3~) is determined 
from eq. OTT l. 

The differential optical depth (|33] | is calculated by separating the contributions of the two resonant populations (3^ and f3~ 
of particles, Afi-es = Afjgs + Af"^. Next we generate a uniformly distributed random number in the range [0,1]. If < 
( 1 - exp [Afres] ), th e photon scatters, in which case the new direction angles and frequency are calculated according to the method 
described in § 13.51 In case both Af""" and Af~^ are non-zero, the value of (3 used for scattering is chosen by generating a new 
random number 5Ri : if 5Ri < Af^^^j,/ Afi-es, the photon scatters with /J"*", otherwise it scatters with f3~. 



3.5. Choice of New Photon Direction after Scattering 

The direction cosine of the photon is chosen randomly in the rest frame^ of the scattering particle, with a weighting appropriate 
to the particular scattering process. We use a similar approach to calculate the surface reflection of an 0-mode photon. 
The differential cross section for resonant cyclotron scattering into both polarization modes is proportional to 

da. 



dfj,'da 

' Througout this subsection, the prime labels the photon after scattering. 



cxl+M'^ (40) 
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(e.g. IWang et aLlll988h . Here fi' = B - k is the direction cosine of the outgoing photon in the rest frame of the scattering particle, 
and a is its azimuthal angle about the local magnetic field direction. The cumulative probability of scattering into an angle less 
than or equal to fi' is then 

1 



' + 3^'+4), 



(41) 



We set P equal to a uniformly distributed random number in the range [0, 1], and then solve for the value of /i'. The direction 
cosine /.i of the photon in the rest frame of the star is obtained from eq. ( fTTI ). since the resonant value of (3 is known. 

The polarization of the outgoing photon is then determined by drawing a second random number in the range [0,1]: the 
outgoing photon is in the E-mode if this number is smaller than (1 +/i'^)"', and is in the O-mode otherwise. 

The azimuthal angle a is uniformly distributed between and 2tt; we take the zero of a to coincide with the B-z plane (see 
Fig. 13.5b . The angle IjLb=B- zis known from the position and magnetic field geometry. We then find the magnetic polar angle /^k 
of the outgoing photon using the cosine theorem for spherical triangles. 



/ik = MB/i+ V'^^~^b)(1~M )cOSQ;. 



(42) 



To find the azimuthal angle 0k about the magnetic axis of the star, we solve for the cartesian components A:^, ky, and fcz, using the 
system of equations 

fcz = Aik (43) 
k■B = ^i (44) 



^.(Bxz)= W(l-/i|)(l-^2)sina. 



(45) 



This gives 0k = arctan(fey/fex). 



3.6. Statistical Treatment of Recorded Data 

The state of each escaping photon is defined by its frequency u; and direction cosine /ik relative to the magnetic axis. (The 
magnetosphere and photon source are axially symmetric, and so we ignore the azimuthal angle 0k.) 

Care must be taken in choosing the frequency distribution of the seed photons. It is straightforward to draw directly from a 
blackbody number distribution 



5 log 



dN 



(Win/Wbb) 



1' 



(46) 



dlog(uJi„/u)bb) exp(win/a;bb)- 

where Win/wbb is the input frequency in units of Wbb = k^^b/fi. However, it is difficult to obtain good sampling of this distribution 
at very low and very high frequencies, because it drops off rapidly at ^ a>bb and lo ^ LO\,b- 

We therefore adopt an alternative procedure: we use the same input frequency for all the photons, and store the output frequen- 
cies and direction angles in a normalized response function 

dN 

(47) 



^1 



(log 












. ^in . 





1 



A'phot dlJ,i^dlog(uJout/uJin) ' 

Here Wout/win is the frequency shift relative to the input frequency, and A^phot is the total number of photons, 
histogram in frequency and direction is then 



//log 



Wbb 



,Mk = / /J log 



/ik B log 



d\og 



i^bb 



The outgoing 



(48) 



which is a convolution in loguj. The underlying assumption in this procedure is that the outgoing distribution in /ik is indepen- 
dent of frequency, which certainly holds for the self-similar field geometry, but not for more complicated field geometries with 
multipolar structure. The drawback of this approach is that by choosing a single input frequency, one implicitly fixes the the 
initial scattering radius (which depends on frequency) relative to the stellar radius. Nevertheless, since most photons are emitted 
at frequencies near the black body peak, excellent accuracy is obtained by choosing the peak frequency 



:LJpeak = 2.82144 



^B^bb 



(49) 



as the input frequency. In Appendix lAl we describe validation runs that confirm that this procedure produces results very close 
to those obtained by drawing frequencies directly from the black body distribution. 

The response function R is normalized to unity (modulo the small effect of photon re-absorption at the surface), because the 
total number of photons is conserved by scattering. The upscattering of photons in frequency therefore implies a transfer of 
kinetic energy from the charges to the photons, which u ltimately is drawn from the energy stored in the external toroidal field 
(frhompson et al. 2002; Beloborodov & Thompso nl2006h . 

We now show how pulse profiles and line-of-sight dependent spectra may be obtained from the function H. The energy flux at 
a fixed magnetic polar angle and in a fixed band tJmin < w < t^max is obtained by integrating H over frequency. 



dF 



}iLuH(Lu,fi]J dloguj. 



(50) 
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To simulate a pulse profile, we pick the angle 9n between the rotation axis and the magnetic axis, and the angle ^los between 
the rotation axis and the light of sight. In a coordinate system that is frozen into the star and aligned with its magnetic axis, the 
rotation axis rotates in azimuth, 

^n = -nt. (51) 

To obtain an accurate pulse profile with minimal noise, we integrate the energy flux over a beam of a finite angular width 
^beam = COS"' /ibeam that is Centered on the direction O to the fiducial observer, 

1 /"^^ /" 1 //jr 

P^=^ ZJ-Qif-' O- Mbeam) dfi'd^' . (52) 

2vr Jo J_i rf^k 

Here, 8 is the step function. The band energy flux so constructed depends only on the rotational phase (pQ and the orientation 
angles Oq, 9\os- The angular integral ( |52] l is performed most easily in an inertial coordinate system that is aligned with the rotation 
axis, and in which O = (^^los,0). Each point r' = (cos"' ii',cj)') in the integral must be related to the magnetic coordinate system, 

= + - - A^'^) cos (0n - 0'), (53) 

and to the observer's direction, 

r'd = Mios/i' + ^(1 - mL)(1 - m") cos 0' . (54) 

In these expressions, 

fiio, = cos(6'ios); fJ.n = cos(6in). (55) 
The rotationally-averaged spectrum is obtained by removing the frequency integral and then integrating over rotational phase, 

^^^ = 7r~3 I I I fi^H(uj,^iy)Q(r' ■d-^be^Lm)dfi'd(l)'d(t)n- (56) 
(47r)2 Jq Jq J_i 

4. RESULTS 

Our main goal in this paper is to determine the basic types of non-thermal spectra and pulse profiles that result from the transfer 
of X-ray photons through the scattering corona of a magnetar To this end, we consider only the simplest distribution of seed 
photons - a single temperature black body - and the simplest example of a twisted magnetosphere, the self-similar sequence 
Q with a magnetic pitch angle that is independent of radius. We only make qualitative comparisons of our calculations to the 
observed population of SGRs and AXPs. For example, features in individual pulse profiles which are not reproduced by our 
calculations could plausibly indicate the presence of hotspots on the stellar surface, or some non-axisymmetry in the current 
flowing out to ~ 1 00 km f rom the star We are able to reproduce the 1-10 keV pulsed fractions that are seen in all but one source 
(IE 1048.1-5957: lTiengo et al. 200 5). 

We do explore in some detail how the X-ray spectra and pulse profiles depend on the momentum distribution f{p) of the 
magnetospheric particles. We consider the two basic types of distribution function itemized in 32.31 the relativistic Boltzmann 
distribution (eq. |fT9l ); and a broad, power-law distribution in momentum (eq. ||2TI ). We also perform some test calculations 
(AppendixlAll with a top-hat velocity distribution ( |36] |. We explore a range of mean drift speeds and power-law indices, and also 
vary the magnetospheric twist angle A0n-s (eq. H]) and the orientation angles /i^ and /iios (eq. 1551 ). The scattering charges are 
taken to have the charge/mass ratio of the electron, with the exception of the simulations of ion scattering presented in § 14.51 The 
results depend weakly on whether the photons are injected in the E-mode or O-mode (§ 14. 6b . 

4.1. A Simple Model of Black Body + High-Energy Power-law 

The X-ray sp ectra of the AXPs can be fi tted by a superposition of a blackbody and a power-law, with a photon index in the 
range^ T = 2—4 (I Woods & Thompsonl2006h . Although a single power-law typically provides a good fit to persistent SGR spectra, 
the Galactic SGR sources are more distant than the AXPs and their spectra are more heavily absorbed. A blackbody component 
is present during transient X-ray afterglows that are observe d following bur sts of intermediate energy (e.g. Ibrahim et al. 2001), 
and in the persistent X-ray emission of SGR 1900+14 (Esp osito et al.ll20()6l) . 

In AXP spectra, the thermal and non-thermal components usually have a similar normalization at an energy of 1 keV. In other 
words, the power-law compone nt pivots about th e black body peak, a behavior that was observed directly during a ~ 1-day 
X-ray outburst of IE 2259+586 (Woods et £01200 4). Within that outburst, which produced more than 100 SGR-like bursts of a 
sub-second duration, the high energy photon index was observed to decay from F = 2 to F = 3.5. 

This behavior suggests that the high-energy extension of the black body peak in the X-ray spectrum is created by upscattering 
the thermal photons in a corona above the surface of the neutron star. The model examined in this paper has this property. It also 
has the advantage that no fine tuning is required for the optical depth of the scattering charges, or their mean energy. The optical 
depth is close to unity if the magnetosphere (or, more precisely, the part of the magnetosphere 100 km distant from the magnetar 
surface) is twisted. 

The standard blackbody + power law fit to magnetar X-ray spectra fails to account for an important feature of our model: the 
high-energy power law component of the spectrum is not expected to continue below the black body peak. This is not just an 

* The photon index V is defined by dN /dui oc uj"^ , where dN/duj is the number of photons emitted per unit frequency. 
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academic issue: in sources with relatively low extinctions, the power-law component dominates below the therm al peak, and 
significandy biases the fit to the black body temperature (e.g. the AXP XTE J1810-197: lHalpern & Gottheljl2005h . 

This leads us to propose a simple generalization of the standard blackbody + power law model. This model assumes that 
a fraction /unscatt of the seed photons leaves the star without scattering, and the rest are upscattered following a power-law 
distribution (Appendix A). The angle averaged response function, eqn. ( |47] |. can be written as 



^ log 



Wout 
Win 



= /unscatt S log 



Wout 

Win 



+/p.e(^^-i 

' Win 



Wout 



r-1 



where O is the step function. Using eqns. ( |46] l and ( |48] l, we get 



Wbb 



+ /pl 



B 



Wbb 



UJ 
UJ 



r-1 



duj' 



(57) 



(58) 



where B(ujluj\jo) is the Planck function, eqn. ( |46] |. Even though this behavior of the response function is seen only in simulations 
with a monoenergetic distribution function, eqn. ( |36] |. this generic fit worked relatively well for most of our results 



4.2. Photon Index 

Figures|4l|6]show the energy spectra (ujF^) that result from choosing particle distribution functions of types I-III, respectively. 
The Planck function is plotted for comparison, and is normalized to have a unit area. A high-energy excess above the black body 
is present, and in most cases is seen to have an approximately power-law shape. 

For each output spectrum, the model described in 34.21 was fit to the data using a Levenberg-Marquardt fitting routine 
(iPress et al.i .l992) in the frequency range -2 < log(a;/a'^™) < 1.7. This fit involved the photon index F, a thermal frequency 
uj^l (which can differ from the input thermal frequency Ci^" = k^Tb^/Ti by up to ^ 10%), and the normalizations of the thermal 
(unscattered) and non-thermal (scattered) terms in eq. ( 15 Sl . These parameters are collected in Tables [T][3] (we only tabulate the 
relative normalization /pi //unscatt)- In some cases, a better fit was obtained by setting /unscatt = in eqn. ( l58b . in which case /pi is 
an overall normalization. In general, most of the fits deviate by about 10% for \og{uj / uj^^) < -1, are accurate to a few percent 
in the range -1 < \og(uj/uj^^) < 1, and deviate by a varying amount (1-100%) at higher frequencies. The label "poor fit" was 
added in Tables 1 and 3 when deviations from the model in the higher-frequency range exceeded ~ 50%. This is usually the case 
when the resonant optical depth is small. 

The hardness of the high-energy spectrum correlates with twist angle, to which the optical depth is directly proportional, and 
with drift velocity of the charge carriers. The dependence of the optical depth on the latter parameter is non-linear, often coupled 
with geometric effects, which results in a stronger dependence of the photon index on the velocity distribution parameters than 
on the twist angle. A mildly relativistic Boltzmann particle distribution function (types I-ll) yields photon indices from F = 2 (for 
Po = 0.8) to F => 7 (for f3o = 0.2 and asymmetric/symmetric particle flows along B). Reducing the twist from A0n-s = 1 results 
in softer spectra. The choice of a broad relativistic distribution (111) produces approximately flat energy spectra above the black 
body peak for particle index a = -l and 7niax = 3, and for a = -2 and a broader range of 7niax- Even harder spectra can be obtained 
by increasing the mean drift speed of the charges, but significant deviations from a power law appear in the high-energy tail. 

A more distinct high-energy component of the spectrum is created when the charges have a broad, relativistic distribution of 
momenta (Fig.|6]l. The energy spectrum dips above the thermal peak and then rises as a power-law up to a maximum frequency 
7iTiaxWpeak- The photou index above the dip can be related dkectly to the index a of the particle distribution by assuming that 
each high-energy photon has scattered only once. The number flux of photons upscattered to a frequency 7^Wpeak is jF^ oc 7/ cx 



1 +r> 
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(l+a)/2 



The high-energy photon index is therefore 



r = 



l-a 



(59) 



A flat particle number distribution (a = -1) results in a rising energy spectrum, wi^j « lo; whereas a flat energy distribution 
(a = -2) results in an energy spectrum that is somewhat softer, ujF^^ oc o;'/^. 

Geometric effects also strongly influence the spectrum that the observer sees. The highly anisotropic character of the emission 
pattern can be appreciated by comparing the spectra of the aligned rotator (/ijj = 1: Figs |4};,e; |5};,e; and Fig.|6^) with those of 
the orthogonal rotator (/ifj = 0: Figs |4jl,f; |5}l,f; and Fig.|6h). Changes in the observer's orientation (the value of /iios) manifest 
themselves through variations in the relative flux of the non-thermal component, and also through changes of the order of AF ^ 1 
in the photon index. Except for the nearly-aligned rotator, the differences in spectral hardness resulting from a change in the 
observer's orientation are smaller than those which divide the spectrum of a quiescent AXP from the spectrum of an SGR. 

4.3. Pulse Profiles and Pulsed Fractions 

The persistent emission of the SGRs and AXPs is pulsed. The pulse profiles generally have one or two main pulses in the 1-10 
keV band, but range from simple sinusoids to more complex, multiple-peaked profiles. The pulsed fraction, defined as 

^max ^min 



PF = 



(60) 



This implies that /unscaii is a fit parameter which may differ from the actual fraction of unscattered photons. 
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TABLE 1 

Spectral Fit Parameters and Pulsed Fractions for Type I Velocity 
Distribution (Boltzmann, unidirectional), A/3 = 0.1 
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" Entri es wi th no data mean that setting /unscait = in eqn )58t provided the best fit.'' Poor fit. See 
section l4^ for details on the accuracy of the fitting procedure. 



TABLE 2 

Spectral Fit Parameters and Pulsed Fractions for Type II Velocity 
Distribution (Boltzmann, bi-directional), A/3 = 0.1 
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" Entri es wi th no data mean that setting /unscatt = in eqn )58t provided the best fit." Poor fit. See 
section |4!2l for details on the accuracy of the fitting procedure. 
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Fig. 4. — Energy spectra obtained with a unidirectional Boltzmann particle distribution (type I; eq. 1191 ). The unit of frequency is o^tb = kgTt,b/n, where the 
black body temperature Tbb is left arbitrary. The dotted line denotes a Planck function of temperature 7i,b> with a unit normalization below the curve. Curves in 
panel (a) correspond to different velocity spreads /3o = {0.2 (short dash), 0.3 (dot-short dash), QA (dot-dot-dash), 0.5 (solid), 0.6 (dash-dash-dot), 0.7 (long dash), 
0.8 (dot-long dash) }. All curves assume that A0n_s = 1 and that the magnetic axis and line of sight are both orthogonal to the rotation axis (/iQ = /tjoj = 0). 
Curves in panel (b) coirespond to different twist angles, A(^n_s = {0.1 (short dash), 0.3 (dot-dash), 0.7 (solid), 1 (dot-dot-dash), 1.3 (dash-dash-dot) }. Velocity 
spread /3o = 0.75 is assumed and other parameters are as in (a). Changing the line of sight direction results in spectra shown in panels (c)-(f): curves correspond 
to Mlo,s = { 1 (dash-dash-dot), 1 /\/2 (dot-dash), (solid), —1 / \/2 (dot-dot-dash), -1 (dash) }. Panels (c) and (e) correspond to an aligned rotator (/in = 1), and 
panels (d) and (f) to an orthogonal rotator Panels (c) and (d) assume a velocity spread /3q = 0.6, and (e) and (f) assume /3o = 0.75. Other parameters are the same 
as in panel (a). Photon indices obtained by fitting a blackbody plus power law function in the frequency range log(aj/a;(,b) = 0.7- 1.4 are listed in Table[T] 

where /^min. max are the minimum and m aximum fluxes measured within one rotation period, ranges from a few percent to as high 
as ^ 70% ( Woods & Thompsonl2006l) . A complex pulse profile is seen more frequently in the persistent emission of the SGRs. 

The pulse profiles resulting from particle velocity distributions of types 1-111 are displayed in Figs.|71|9] (These figures conre- 
spond to the spectra in Figs. |4l|6]) The profiles are plotted in five frequency bands, which are taken to be 
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(61) 



where Wpeak is defined in eq. ( |49] l. Each pulse profile is normalized by the pulse-averaged flux in that band: averaging over the 
rotation gives unit flux. The pulsed fractions are listed band by band in Tables [T]|3j the values are generally accurate to 1%. 
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The pulse profiles have a wide variety of morphologies, and sometimes vary significantly with frequency. Several properties 
can be understood from the fact that the direction of a scattered photon is collimated more tightly about the direction of the 
magnetic field as the speed of the scattering charge increases. The differential cross section for resonant scattering is proportional 
to (1 + /i'^), where ^i' is the direction cosine in the rest frame of the charge after scattering (e.g. Meszaros 1992). The probability 
of scattering at an angle /i > in the stellar frame is Pip. > 0) = (4+3/3+/3-')/8 [see § 13.51 . For (3 = 0.5, there is a 70% probability 
of photons being scattered in the same direction as the charge is moving. 

Several trends are apparent in Figs.|7]|9] The pulsed fraction increases with particle drift speed, as does the asymmetry between 
the main pulse and the sub-pulse 180 degrees away in phase (Fig.|7^). The main pulse is emitted toward the south magnetic pole, 
and the sub-pulse toward the north magnetic pole. Although the scattering particles flow out of the north pole and into the south 
pole, the scattering depth is largest near the magnetic equator, where the particles flow southward. The asymmetry between main 
pulse and the sub-pulse also grows with increasing twist, because of the growing asymmetry in the photon flux near the magnetic 
equator due to the increased optical depth there (Fig.lTJ?). At the same time the main pulse divides into two sub-pulses, because 
the density of charge carriers vanishes along the magnetic axis. The sub-pulse is generally more distinct in the bands near the 
peak of the input thermal spectrum, because these photons are less scattered that those at higher energies (Figs.|7};,e). Finally, the 
pulsed fraction is generally larger for nearly orthogonal rotators (Figs. |7}l,f). 

In spite of this diversity, the axisymmetry of the magnetic field translates into a reflection symmetry of the pulse profile about 
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Fig. 6. — Energy spectra obtained with the broad momentum distribution, eq. 12 It . Normalization and frequency units are the same as Fig.|4] Panel (a) shows 
results obtained with different /3n,i„, for "fmux = 3, a = —1, A0n-s = ~U and an orthogonal configuration, fiQ = ii[„s = 0. Curves correspond to /3n,i„ = {0.1 (short 
dash), 0.3 (solid), 0.7 (dot-dash), 0.9 (long dash) }. Panel (b) shows the result of varying the twist angle A</<n_s> for /^min = 0.3, 7max = 30, a = -l, and the same 
orthogonal configuration as (a). Curves coiTespond to A0n_s = {01 (short dash), 0.3 (dot-dash), 0.7 (solid), 1 (dot-dot-dash), 1.3 (long dash) }. Panel (c) shows 
the result of varying 7max, for /3,„in = 0.3 and identical parameters as (a). Curves correspond to 7max = {1.1 (short dash), 1.96 (solid), 3.48 (dot-dot-dash), 6.19 
(dash-dash-dotted), 11 (dot-long dash), 19.6 (dot-short dash), 34.8 (long dash) }. This calculation is repeated in panel (e) for a softer particle distiibution with 
index a = —2. Panel (d) shows results obtained by varying the exponent in the distribution function a, for 7niax = 30 and other parameters the same as in (c). 
Curves correspond to a = {—4 (short dash), —3 (dot-short dash), —2 (solid), —1 (dot-dot-dash), (dot-long dash)}. This calculation is repeated in panel (f) for 
7max = 3. Finally, panels (g) and (h) show the aligned and orthogonal rotator configurations, respectively, analogous to those in Fig.|4](c)-(f), but with simulation 
parameters /3,^in = 0.3, 7mM = 30, a = -1, and AipM-s = 1- 
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TABLE 3 

Spectral Fit Parameters and Pulsed Fractions for Type III Velocity Distribution 
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^ Entries with no data mean that setting /unscatt = in eqn )58t provided the best fit.'' Poor fit. See section I431 for 
details on the accuracy of the fitting procedure. 



phase = TT. The absence of such a mirror symmetry is a direct signature of some non-axisymmetric structure in the surface 
heat flux or in the magnetospheric currents, which is not accounted for in the present study. 

The pulse profile also provides information about the composition of the magnetospheric corona. Some of the trends in Fig.|2] 
are preserved when the velocity distribution is of type 11 (symmetric under /3 — > -/?), but some differences emerge. Figure [8] 
displays the resulting pulse profiles; all other parameters are held constant from Fig. |7] When the magnetic axis is nearly 
orthogonal to the rotation axis (as in panels a,b,c and e), the observer must see two comparable pulses, because the charges 
stream in both directions along the magnetic field. For an orthogonal rotator, as displayed here, the light curve between rotational 
phase 180° and 360° is, in fact, a mirror image of the light curve between 0° and 180°; but this orientation is a special case. More 
generally is possible for the observer to see a single pulse, as demonstrated in Fig.|8}l,f in the case where /io = 1/V2- 

A single pulse can divide into two sub-pulses for the same reason as in Fig.|7] the current density and scattering rate decreases 
toward the magnetic axis. One interesting feature of Fig.[8]is that the positions of the sub-pulses can alternate by tt/2 in phase 
between different bands: this reflects the number of scatterings that a photon in each band has typically experienced. 

It should be emphasized that, aside from the number of major sub-pulses, which depends directly on the composition of the 
magnetospheric plasma, most properties of the pulse profiles depend even more strongly on orientation. A good example is 
provided by Fig. |7f, which shows how a single-peaked profile with subpulses (dotted line) can be transformed into a profile 
without any subpulse (dashed line), and then to a double -peaked profile (solid line) merely by changing fi^ and /iios- 

The calculated pulsed fractions P range from 0% to nearly 100%. The magnitude of P and its dependence on the simulation 
parameters can be quite different in the low-energy and high-energy bands. In general, the pulsed fraction correlates positively 
with twist angle as long as A(/)n-s ^ 1- The largest pulsed fractions are obtained for the broad velocity distribution (eq. ET\ ). as 
is shown in Fig.|9] 

It should be kept in mind that the heat flux is expected to vary strongly over the surface if the elastic deformations of the 
crust are loca lized in a small are a : the small radiative area s of SGR and AXP afterglows (about 1 percent of the neutron star 
surface area: llbrahim et al.ll200l1; IWoods & ThompsOTill200 6) provide evidence for this behavior Our assumption of a single- 
temperature blackbody therefore involves a considerable simplification of reality, at least in the more transient magnetar sources. 
The combined effects of anisotropic emission and magnetospheric scattering are not examined here. 

4.4. Spectral Features 

iLvutikov & Gavriil (1200 6') suggest that spectral features in the seed photon frequency distribution may be smeared out as the 
result of multiple scattering in a magnetar magnetosphere. Here we test this effect with our three-dimensional Monte Carlo 
model. 

Figure [To] shows energy spectra obtained with a narrow velocity distribution (/3 = 0.5, A/3 = 0.1) and large twist (A(/)n_s = 1). 
Instead of convolving the response function with a plain black body spectrum, we now add an emission line to the seed photon 
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Fig. 7. — Pulse profiles obtained with a unidirectional Boltzmann particle distribution (type I; eq. 1191 ). All curves are normalized so that the av erage amplitude 
is unity. Simulation parameters of (a), (b), (c)-(d), (e)-(f) are the same as in the corresponding panels in Fig.|4] Energy bands are defined in eq. )6U . Curves in 
panel (a) show the band V profile for different velocity spreads: /3o = {0.4 (dash), 0.5 (dot-dash), 0.6 (dotted), 0.7 (dash-dash-dot), 0.8 (soUd)}. Band V pulse 
profiles in panel (b) correspond to different twist angles: A0n-s = {0.1 (dash), 0.3 (dot-dash), 0.7 (dotted), 1 (dash-dash-dotted), 1.3 (solid)}. Curves in panels 
(c) and (e) are for different bands: I (dotted), II (dot-dash). III (solid), and IV (dash), with /3o = 0.6 and 0.75, respectively. Band V pulse profiles in panels (d) 
and (f) correspond to different orientations of the magnetic axis and the line of sight with respect to the spin axis: (/xn , /.tios) = {(1/ V2, 1 / \/2) (solid), (1/ \/2, 0) 
(dash), (0,0) (dotted)}, with /3o = 0.6 and 0.75, respectively. Pulsed fractions in bands I-V for all configurations are listed in Table[T] 

distribution. The line is centered at = S^BTbh/^, with half-width Q.\LO\,b, and total flux 1% of the integrated blackbody flux, 
as done by Lvutikov & Gavriil (2006). The spectra observed from both an aligned rotator and an orthogonal rotator are plotted. 
Even though the line strength is reduced by a factor 2 in some orientations, the line is not smeared out. The basic reasons 
are that i) the optical depth is never very large along any direction; and ii) a substantial fraction (of the order of one half) of 
the emitted black body photons must escape the magnetosphere directly without scattering, even if the magnetospheric twist is 
strong. The greatest amount of smearing would be expected in a geometry where a significant optical depth is maintained along 
the line of sight over an entire rotation of the star. Nonetheless, a significant line feature remains even in the case of a nearly 
aligned rotator, with the line of sight nearly orthogonal to the magnetic and rotation axes (the lowest curve in Fig.fTOk). 



4.5. Ion Scattering 

We have performed some preliminary calculations of the effects of ion cyclotron scattering on the output spectrum. We consider 
the case of protons (with a charge/mass ratio 1 /1836 that of the electron). A polar magnetic field of 1 x 10'^ G is chosen, which 
allows proton cyclotron scattering over an order-of-magnitude range in frequency above the black body peak. 
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Fig. 8. — Same as Fig.|4] but for a velocity distribution function that is bi-directional, corresponding to a pair plasma (t™e II; eq. II9I ). Simulation parameters 
of (a), (b), (c)-(d), (e)-(f) are the same as those in the corresponding panels of Fig.|5] Pulsed fractions are listed in Table|2] 

Energy spectra and pulse profiles are shown in Fig.[TT| The results for protons are plotted side by side with those for electrons, 
using the same position-independent velocity distribution (eqs. [T91I2TI) . and varying the angular distributions of the seed photons. 
Solid lines and dot-dashed lines show results for protons assuming radial and isotropic seed photon emission, respectively, with 
solid and short-dashed curves containing correponding results for electrons. The first two sets of figures assume a narrow top 
hat velocity distribution. The case of a bi-directional proton velocity distribution can be related to a situation in which the X-ray 
flux at the proton cyclotron line e xceeds the Eddington fl ux, and a cloud of ions and electrons is blown off the stellar surface and 
suspended in the magnetosphere ([Thompson et al.ll2002h . The final set of figures assume a broad particle momentum distribution 
(type III; eq. ED). 

Multiple ion scattering generates a high-frequency power law tail to the input thermal spectrum. The spectrum is however 
softer than that produced by electron scattering: photons that are backscattered by ions to the star can be absorbed at the stellar 
surface (with unit probability in the E-mode, and with probability ^ | in the O-mode). As expected, the high-energy spectral 
tail is cut off sharply at a frequency (rbb/^Ns)*^^''^Wpeak ^ lOwpeak- In this calculation, the cutoff frequency for electron scattering 
is nip/me ^ 10^ times higher We have not displayed this cutoff because the calculation does not include the effects of electron 
recoil (which becomes important at 50-100 keV), or relativistic corrections to the scattering cross section at the first Landau 
resonance. Recoil effects can be entirely neglected for ion scattering. 

The effects of ion scattering could in principle be discerned during a very luminous phase, e.g. the afterglow following an 
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Fig. 9. — Pulse profiles obtained with the broad velocity distribution function, eq. 12 IK with same normalization as Fig.|2] Simulations parameters are as shown 
on each panel. Panel (a) shows profiles in band V obtained with different minimum velocities, with curves corresponding to = {0.1 (dash), 0.3 (dotted), 
0.7 (dot-dash), 0.9 (sohd)}. Panel (b) shows the result of varying the maximum Lorentz factor, with curves corresponding to 7n,ax = {11 (dash-dash-dot), 1.96 
(dot-dot-dash), 3.48 (dash-dot), 6.19 (long dash), 11 (dotted), 19.6 (short dash), 34.8 (solid)}. Panel (c) shows results for different twist angles, with curves 
corresponding to A(/iim_s = {0.1 (dash-dash-dot), 0.3 (dash-dot), 0.7 (dotted), 1 (long dash), 1.3 (solid)}. Panel (d) shows the result of varying the exponent of 
the distribution function, with curves corresponding to q = {—4 (dot-dash), -3 (short dash), —2 (solid), —1 (long dash), (dotted)}. Panel (e) shows lighcurves 
in bands I (dash), II (dotted). III (dot-dash), and IV (solid), for /3,ni„ = 0.3, 7niax = 3, « = — 1, A0n-s = 1> and = ii\as = 0. Panel (f) shows profiles in band V 
obtained from the same simulation as in (e), but vaiying the orientation as (/in,/iios) = {(l/\/2, l/\/2) (solid line), (l/\/2,0) (dashed line), (0,0) (dotted line)}. 
Pulsed fractions are shown in Table|3] 

intermediate-energy SGR burst. For example, the afterglow of the 29 August 1998 burst of SGR 1900-1-14 showed a marked 
spectral hardening in the 2-5 and 5-10 keV bands, but not in the 10-20 keV band (I brahim et al. 2001.) . 

Results for proton scattering are sensitive to the angular distribution of seed photons, given their smaller resonant radius 
(Fig. fTTI) . Their spectra become harder when the angular distribution is changed from radial to isotropic. In contrast, spectra 
obtained with electron scattering are almost independent of the input angular distribution, with results for the radial and isotropic 
case being hardly distinguishable. 

The pulses produced by ion and electron scattering have a similar morphology (Fig. fTTI) . although electron scattering yields 
smaller pulsed fractions. Again, results for protons are somewhat sensitive to the angular distribution of seed photons, while 
the same parameter does not influence pulse profiles obtained with electron scattering. A proper account of ion scattering must 
include the effects of light bending of the photon trajectories 
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Fig. 10. — Effect of magnetospheric scattering on an emission line. Input spectrum (solid curves) has a line centered at a; = 5k^Ti,b/H, with half-width 
OAkBTbb/h and total energy equal to 1% of the blackbody. Output spectra are obtained by convolving the input with the magnetospheric response function for 
different orientations. Left: Aligned rotator (/in = 1), with /iios = 1 (dash), 1 / Vl (dotted), (dot-dash), -1 /\/2 (dot-dot-dash), and -1 (dash-dash-dotted). Right: 
Orthogonal rotator {/iq = 0), with ^los = 1 (dotted), 1 /\/2 (dash), and (dot-dash). Velocity distribution is mono-energetic and undirectional (eq. 1361 ). with a 
mean drift speed /3 = 0.5 and a total width A/3 = 0.1. The twist angle is A4>n-S = 1- 

4.6. Seed Photon Polarization 

There are two possible sources of seed X-ray p hotons: internal heating by a decaying magnetic field (e.g. iThornpson & Duncanl 
19961: iHeyl & Kulk:arnil[l998l: lArras et aP 120 04'). and the bombardment of the surface by energetic charges dThom pson et al.l 
2002h . In the first case, the seed photons are emitted in the E-mode if the magnetic field is stronger than '-^ 10^^ G (Ho & LaH 
2004). In the second case, the emission is in the O-mode if the charges are stopped by Coulomb collisions in a surface layer that 
is optically thick to free-free absorption. 

It is therefore interesting to explore whether these two types of seed photons can be distinguished using the output X-ray 
spectrum or pulse profiles. The scattering cross sections of the two modes at the cyclotron resonance are comparable, because 
the electromagnetic eigenmodes are determined by vacuum polarization effects and are linearly polarized even in the core of the 
cyclotron resonance (©. For this reason, one expects similar pulse profiles and non-thermal spectra for both polarizations of 
the seed photons. Nonetheless, the scattering cross section of an O-mode photon is suppressed with respect to that of an E-mode 
photon by a factor (/i')^ (eqs. fTSl . ITSl : see also the discussion in Appendix lAli. One therefore expects the spectra to be slightly 
softer and the pulsed fractions to be slightly lower when the seed photons are injected in the O-mode. 

Figure[T2]shows energy spectra (left) and pulse profiles in band V (right) using the same magnetospheric parameters, excepting 
that the seed photons are 100% polarized in the E-mode (solid lines) and the O-mode (dashed lines). In the second case, the non- 
thermal tail does indeed have a slightly lower amplitude compared with the black body peak, and the pulsed fraction is smaller 
Nevertheless, the changes in the emission pattern resulting from the polarization of the seed photons are much subtler than those 
produced by changes in the orientation of the star ( ^4.21 14.3b . offering little chance of identifying the seed photon emission 
mechanism except through direct polarization measurements. 

5. DISCUSSION 

We have developed a Monte Carlo co de to study the reprocess ing of thermal X-rays by resonant cyclotron scattering in neutron 
star magnetospheres using the model ofi Thompson et all(l2002h . The code is fully three dimensional, and allows for an arbitrary 
velocity distribution of the scattering particles and an arbitrary magnetic field geometry. Angle-dependent X-ray spectra and 
pulse profiles can be calculated for an arbitrary orientation between the magnetic axis and the spin axis, and between the spin 
axis and the line of sight to the star 

Our aim here has been to show which types of X-ray spectra and pulse profiles follow generically from the model, and which 
may point to different emission mechanisms. We can reproduce most of the properties of the 1-10 keV emission of the AXPs, but 
not of the SGRs, by assuming a broad and mildly relativistic velocity distribution, and moderate twist angles (A0n-s ~ 0.3- 1 
radian). In p articular, the SGRs in their most active states have harder X-ray spectra than can be reproduced by the model 
(I Woods et al.ll2006i) . Our results demonstrate that changes in the strength of the twist and in the dynamics of the charged particles 
(which are coupled strongly to the radiation field) can lead to substantial changes in flux, hardness, pulse shape, and pulsed 
fraction. 

5.1. Relative Strength of the Thermal and Non-thermal Components of the X-ray Spectrum 

The X-ray spectra of most magnetars display a prominent thermal compone nt. The cor responding blackbody area is less than, 
or sometimes comparable to, the expected surface area of a neutron star (e.g. IWoods & Thompson 2006). In sources where the 
blackbody component is dominant at ^ 1-3 keV, the high-energy photon index is typically 2 or softer SGR 1900+14 presents 
an example of an SG R where a thermal component and a hard power-law component of the spectrum are simultaneously present 
dEsposito et al.ll2006l) . 
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Fig. 1 1 . — Energy spectra (left) and band V pulse profiles (right) obtained for Bpoie = lO'^ G and a chai'ge/mass ratio equal to that of the proton and the electron, 
for radial (long-dash: proton, solid: electron) and isotropic (dot-das hed: proton, short-dashed: electron) angular distribution of the seed photons. Panels (a) and 
(b) are for a mono-energetic, unidirectional particle distribution (eq. 1361 ) with mean drift /3 = 0.7 and A/3 = 0.1. In panels (c) and (d), the velocity distribution is 
mono-energetic and bi-directional, as is appropriate to the case where the ions ai'e driven oft" the s urfa ce of the magnetar by radiation pressure. Other parameters 
are the same. In panels (e) and (f), the velocity distribution is broad and relativistic (type III; eq. 1211 ). with = 0.3, 7max = 3 and a = -1. All curves assume 
'^0N-s = 1> and MO = A'los = 0. Curves for electrons are almost undistinguishable. 

Model spectra with a mildly relativistic particle distribution (ksTo ~ 0.5meC^,/3o — 0.75; eq. |[T9l ) can easily reproduce the 
1-10 keV spectra of sources like 4U 0142-1-614, IRXS J170849-4009 and IE 1841-045. Some of our model spectra show 
a noticeable downward break at ~ 30-50 A:Brbb — 10-30 keV (e.g. Figs. 4f, 5f). Such a break may be present in the combined 
XMM and Integral spectra of SGR 1900-1-14 ( Gotz et al. 2006), and its presence is not excluded in other sources such as the AXPs 
4U 0142-1-614 a nd IE 1841-045 when one takes in to account the presence of a separate rising high-energy spectral component 
above 20 keV ( Kuiper et al.ll2006t iGotz et sIMo^ . 

The quiescent SGRs have relatively weak blackbody components t o their spectra, and during periods of extr eme activity a 
pure power-law fit to the spectrum can have a photon index of ~ -1 .6 jMereghetti et alJl2005l : IWoods et al.ll2006 !). (Even harder 
photon indices are obtained from combined blackbody-powerlaw fits to the spectrum, but then the continuation of the 2-10 keV 
fit to higher energies disagrees with the Integral flux measurements.) We could obtain spectra as hard as this by choosing a 
broad power-law momentum distri bution for the m agnetospheric charges (Fig. 6), but the output spectrum then displays a strong 
blackbody component. In the fit of lEsposito et al.l (2006) for SGR 1900-1-14, the amplitude of the black body component is about 
three times that of the power-law component at the black body peak, whereas in the spectra of Fig. 6b, it is 6-8 times larger. 

We expect this conclusion to hold for any model for the high-energy continuum which relies on the upscattering of the thermal 
seed photons. Balancing the energy input to the coronal charges with the loss by upscattering, one finds that r((7^) - 1) '--^ 1 (as 
in corona models based on non-magnetic Compton scattering). When (7^) '/^ ^ 1, it is generally not possible to upscatter a large 
fraction of the thermal photons because the optical depth r of the corona is small. The results obtained with a broad relativistic 
particle distribution function should also be treated with caution: the charges are subject to a strong cyclotron drag force at a 
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Fig. 12. — Energy spectra (left) and band V pulse profile (right) corresponding to 100% pol ariz ation of the seed photons in the E-mode (solid curves) and the 
0-mode (dashed curves). Velocity distribution is mono-energetic and and unidirectional (eq. (36)) with mean drift speed /3 = 0.7 and total width A/3 = 0.1. In 
addition, A(/>n_s = 1, and /iq = /xios = 0. 

radius of ^ 50-100 km where the cyclotron energy is 1-10 keV, and so it is probably not self-consistent to assume that they have 
the same distribution function there as they do in the inner parts of the magnetosphere. 

These considerations point to the conclusion that the hard non-thermal spectral states of the SGRs (with photon indices < 2) 
have a similar origin to the excess high energy emission that is observed from several AXPs above 20 ke V. In the magnetar model, 
this high-energy emission must involve a different emission process, such as bremsstrahlung from a hot transition layer between 
the relativistic corona and the colder neutron star atmospher e; or synchrotr on emission triggered by pair cascades at a distance of 
^ 100 km, where the electron cyclotron energy is ^ 1 keV (iThompson & Beloborodov 2005). 

Indeed, there are reasons to expect that some active magnetars (e.g. the SGRs) can have significantly reduced X-ray emission 
due to deep internal cooling, even when the time-averaged rate of magnetic field decay remains high. Large increases in the 
neutrino emissivity are expec ted when the core neutrons undergo a superfluid transition, which leads to an order-of-magnitude 
drop in blackbody X-ray flux ( Arras et al.l2004 l). The greater intermittency of the magnetic dissipation in a Soft Gamma Repeater 
(as compared with the AXPs) could then be tied to the greater rigidity of the neutron star crust at lower temperatures. Given the 
high absorbing columns Nh ^ 2-6 x 10^^ cm~^ that are inferred from the soft X-ray spectra of the two active SGRs, it would 
then be possible to hide the surface thermal component from detection. (The combination of a low thermal luminosity and 
strong absorption does not remove the inconsistency with the model spectra in Fig. 6c-e: the models then would predict that the 
non-thermal continuum is much dimmer than ^ 10 ergs s"' at a photon energy of a few keV.) 

5.2. Low Energy X-ray Spectra 

Our results show that the rescattering of thermal photons by magnetospheric charges does not lead to any enhanced emission 
at energies below the peak of the seed X-ray distribution. Our new fitting formula, eq. ( fSSl l. takes this low energy behavior of 
resonant scattering into account. 

It should be kept in mind that the softest power-law components in measured spectra are typically needed to fit a low-energy 
excess to a single-temperature black body. In our model most of the emission below the thermal peak comes from unscattered 
photons: our rotationally-averaged spectra have a Rayleigh- Jeans shape at low energies, which lies just below the input blackbody 
distribution. We therefore suggest that a fit to an AXP X-ray spectrum involving a small Nh and a very soft power law component 
is physically inconsistent, and that in such a case the fitting formula dSSl ) is strongly preferred. Of course, the softest AXP spectra 
can also be fitted by a superposition of two or more black body components (e.g. Halpern & Gotthelf 2005). 

The formation of a high-energy tail to the black-body spectrum is closely tied, in this model, to an increased spindown rate (due 
to the flaring out of the poloidal field lines a s they are twisted). The AXP IE 2259+586 provides an example of a magnetar which 
has demonstrated SGR-like burst activity (IWoods et al.ll2004h . but which also has the lowest spindown rate of any magnetar 
candidate, as well as very sof t 1-10 ke V spectrum, and only an upper bound on the 100 keV emission that has been detected 
from several other AXPs ( Kui per et alj |2006). We would therefore expect that the 1-10 keV spectrum of IE 2259+586 is most 
accurately modelled by a variable surface temperature, rather than by a combination of thermal and non-thermal components. 

Indeed, a recent reanalysis of the X-ray spectra of several AXPs h as revealed a broader spectrum than a pure single-temperature 
black body near the thermal peak dPurant & van Kerkwiikll2006ah . This suggests that the thermal X-ray flux is distributed in- 
homogeneously across the surface of the neutron star. This is expected given the inhomogeneous nature of the yielding be- 
havior in the crust of a magnetar, and the a nisotro py of the thermal condu ctivity in the presen ce of a strong magnetic field 
(IThompson & DuncanlfTgMJLvubarskv et al.ir200a iPerez-Azorm et al.ll2006l: lOeppert et al.ll2005h . 

5.3. Pulse Profiles and Pulsed Fractions 

Our calculations generally reproduce X-ray pulse profiles that are characteristic of the AXPs: single or multiple peaks, with 
sinusoidal shapes and sometimes the presence of narrow sub-pulses. The strong angular dependence of the optical depth to 
resonant scattering explains why phase resolved spectra show variations both in their flux level, as well as in hardness (e.g. 
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iTiengo et al.ll2005 l: iGohler et al.ll 20()5h. We find that the X-ray spec trum is generally hardest at pulse maximum, in agreement 
with some measurements (Meregh etti et al.ll2004tlTiengo et alJl2005h . 

The number of sub-pulses depends essentially on the orientation of the spin and magnetic axes of the neutron star, which are 
not constrained by any independent measurements. The pulse profiles of the actively bursting SGRs tend to be more complicated, 
which may reflect presence of high-order multipoles in the current flowing near the stellar surface. As discussed in ^5.11 their very 
hard X-ray spectra are more consistent with a separate emission process, possibly localized close to the surface of the neutron 
star. 

Pulsed fractions as high as ^ 40 - 50% near the black body peak are easily obtained from the anisotropic scattering of ther- 
mal photons in the magnetosphere (Tables 1-3). Significantly higher pulsed fractions are possible in the power-law tail of the 
spectrum, due to the increased anisotropy of the multiply-scattered photons. It should be emphasized that this result is obtained 
assuming a spherically symmetric source of seed photons, i.e. a constant temperature across the surface of the neutr on star Only 
in one object (the AXP IE 1048. 1-5957) is the pulsed fraction significantly higher (70-80%) in the 1-10 keV range (iTiengo et alJ 
l2005h . It should be possible to obtain higher pulsed fractions if the seed photons are emitted at localized hotspots, as suggested 
by the small black bo d y emitting area ([Duran t & van Kerkwiik 2006b). Our calculations are therefore complementary to those 
of iPsaltis et ai] (l2000h . lOzel etd] ( l200ir and lOzell (I2002l) . who studied the pulse profiles and pulsed fractions that result from 
isolated hotspots on rotating neutron stars, and compared their calculations with the observed 1-10 keV pulse profiles of the 
AXPs. 

5.4. Comparison with Previous Work 

Semi-analytical work on resonant cyclotron scattering in the atmospheres of neutron stars goes back to the 1980 s, the orig- 
inal focus being the transfer of radiation through the atm osphere of the star, or through an accretion column ( Yahell fr979t 
IWasserman & SalpeteiifT980l: lNageilfT980l: Im eszaros & Nag el 1985a b). In these calculations, the geometrical thickness of the 
scattering material was generally taken to be small compared with the gradient scale of the magnetic field. The effect of mag- 
netic field inhomogeneity was analyzed by Zheleznyakov (1 996), who found solutions to the one-dimensional radiative transfer 
equation using the Schwarshild-Shuster method. ILyutikov & Gavriill (l2006h generalized this method to include the combined 
effects of the particle motion and magnetic field inhomogeneity. In the limit of small particle dispersion /3c, they found that the 
radiation spectrum that is transferred through an optically thick medium is upscattered in energy by a factor 1+2/3. This model 
was shown to provide a good fit to the 1-10 keV spectrum of IE 1048.1-5937, the AXP source in which the thermal component 
of the spectrum peaks at the highest frequency. However, this non-relativistic and one-dimensional approximation cannot provide 
a realistic description of the more extended high-energy spectra that are observed in some magnetars, including the slope of a 
power-law tail to the input thermal spectrum, the relative amplitude of the power-law and thermal components, or the detailed 
shape and frequency-dependence of the pulse profiles. 

Monte Carlo codes for resonant cyclotron scattering through the surface layers of a neutron star were developed around the 
same time (e.g. Pravdo & Bussard 1981; Wang et al. 1988), and also applied to the problem of line formation in gamma-ray 
bursts (I Wang et al.lll993l) . The last two studies remain the most realistic to date, as they include the effects of thermal broadening, 
vacuum and plasma polarization, as well as electron recoil. All of these effects but the last one are included in our study (plasma 
makes a negligible contribution to the dielectric properties of the magnetosphere at frequencies of 10^^- 10'^ Hz). However, 
these studies are all one-dimensional and focus on the regime of very large optical depth and a weak gradient in the magnetic 
field, and so the results cannot be compared du^ectly with ours. 

5.5. Further Developments 
We now discuss how some of the key simplifications of the model will need to be relaxed in future work. 

1. Anisotropic radiation pressure. The most important simplification in the treatment of the transfer of radiation through 
the magnetosphere is the decoupling of the particle dynamics from the radiation force that is imparted through cyclotron 
scattering. An electric field is induced by the twisted m agnetic field l ines, which lifts ch arged particles off the stellar 
surface, and accelerates them to high speeds ([Thompson et al. 2002, Belo borodov & Thomps on 2006). Where the Lorentz- 
boosted cyclotron energy lies in the range 1-10 keV (or 1-100 keV for an active SGR), the particle moving against the 
direction in which the radiation i s flowing feel a strong drag force, and a strong electric field is required to compensate 
([Thompson & Beloborodovl2005h . The charged particle dynamics and the radiation transfer are therefore strongly coupled. 

2. Particle recoil. The non-inclusion of particle recoil after scattering leads to significant errors in the shape of the output 
spectra above ~ 100 keV. The spectra obtained with broad relativistic particle distributions involve essentially the single 
scattering of photons of a frequency ^ Wpeak = 2.82^B7bb/^- Recoil is important only when 7 > meC^/2/iajpeak and hu ~ 
T^fjJpeak ^10 MeV, where we do not tabulate spectra. It should be emphasized that this remains a technical issue insofar as 
additional emission mechanisms dominate at ^ 30-100 keV (see 35.1b . 

3. Seed photon distribution. In all our simulations we have assumed that the seed photons are emitted radially, from points 
uniformly distributed over the stellar surface, and with a single polarization. There are several reasons for this choice. 
First, the radius of first scattering by electrons (or positrons) is /?res ^ 10/?ns; second, we are not including the effects 
temperature inhomogeneities on the neutron star surface, focusing instead on in the pulsations that result from the transport 
of X-ray photons through a twisted magnetosphere. The effects of light bending must be included as soon as this second 
assumption is relaxed. 



26 Fernandez & Thompson 

4. Polarization. Continuum spectra and pulsed profiles do not provide a clear way of discriminating between thermal emission 
from the surface due to deep cooling (E-mode emission) and due to particle bombardment and stopping by Coulomb 
collisions (O-mode emission). In both cases, the surface temperature will in practice be anisotropic (due to the channeling 
of the cooling heat flux along the magnetic field, and due to variations in the strength of the twist across the surface). In 
the second case, one expects the inten sity of the O-mode photons to be more strongly beamed along the local direction of 
the surface magnetic field (e.g. B asko & Sunyaevlll97 5h. but we have found that pulsed fractions of 40-50% percent are 
achievable solely by radiation transport through the magnetosphere. Polarization measurements of AXPs and SGRs would 
say much about the physics of the beam-heated atmosphere of a magnetar, and the distribution of magnetospheric currents, 
but they require polarimeters that are sensitive at flux levels of ^ 10"'' ergs cm~^ s~' in the 1-10 keV band. 

5. Magnetic field geometry. The magnetospheric model employed by [Thompson et all (120021) is self-similar: the pitch angle 
of the external field lines is independent of distance from the neutron star Clearly different X-ray spectra and pulse 
profiles will result if the current is distributed in a different way throughout the magn etosphere. Observat ions which show 
a delayed increase in the spindown torque following periods of X-ray burst activity (IWoods et al.ll2002h suggest that the 
twist is not distributed homogeneously throughout the magnetosphere in the immediate aftermath of a crustal yielding 
event. Quadrupole and octopole components of the magnetic field have a large effect on the pulse profile only if ions are 
the dominant scattering s pecies, or if the no n-thermal emission is localized close to the neutron star surface (as it is with 
bremsstrahlung emission: [Thompson & Beloborodov.2005i) . 
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APPENDIX 
A. CODE TESTS 
Input Frequency Distribution and General Output 

The output spectrum resulting from the input of a monochromatic line (frequency Uin) is obtained by integrating the response 
function R (eq. p7l ) over magnetic polar angle /ijj. Well sampled spectra and pulse profiles are obtained with 10^ photons per 
run, using bins of size 0.02 in log(ti>)_and 0.01 in /ik. The left two panels of Fig. [T3] show the result for a mono-energetic particle 
distribution function (eq. 1361 ) with f3 = 0.5, A/3 = 0.1, and particle flow in either one or both directions along the magnetic field. 
The two curves in each graph represent the injection of the seed photons in the two polarization eigenmodes. 

The output spectrum has a sharp peak at ujom = ^in- The fraction of the seed photons that escape unscattered is shown in the 
right panel of Fig. [T3] It is higher for the O-mode than for the E-mode, because the O-mode cross section has an additional factor 
of in')^, where iJ,' = k-B is evaluated in the rest frame of the charge (eqs. ifTSl , ifTSl ). The fraction of unscattered photons is 
highest along the two magnetic poles, where the optical depth goes to zero. Here it is independent of the polarization of the seed 
photons, because the two eigenmodes are degenerate when a photon propagates nearly parallel to B. 

The fraction of unscattered O-mode photons shows a third peak, which arises because the cyclotron optical depth vanishes 
where p, = j3. This condition is equivalent to /i' = 0: the photon propagates orthogonally to the magnetic field in the rest frame of 
the charge. (It is satisfied at /i^ — 2.77 for radially moving photons in a dipole magnetic field, where the particle drift is in only 
one direction along the field at a speed (3 = 0.5. The small bump in the E-mode curve at this angle is a numerical artifact.) 

Above the input frequency, the escaping photons have an approximately power law distribution. The power-law tail is present 
because i) the optical depth r through the corona is independent of resonant frequency (eq. fTTI ). and ii) because the particle 
distribution function is assumed to be independent of position. It follows from i) that each photon has a probability 

P(«)= [l-exp(-T)]" (Al) 

of undergoing n scatterings. At the same time, ii) implies that the frequency increment Aw due to backscattering is independent 
of the scattering order. The probability that a photon is upscattered from a frequency Win to w,, = (1 + Acj/(jj)"win is therefore 

/ ^ {ln[l-exp(-T)]/ln[l+Aw/ii;]} 

P(uj„)=(^] . (A2) 

The power-law index is negative here because Ao; > 0. It is of the order of unity when the particle drift is mildly relativistic, 
because Atj/w ^ 1. A more precise analytic cal culation of the spectral slo pe is very difficult in practice, because r and Alj/cj 
depend on p, and /3 in a non-linear manner. See iPozdnvakov et alj (11983b for a general review of the effects of multiple non- 
resonant Compton scattering in X-ray spectra. 

The spectral signature of the first few orders of scattering is apparent in Fig. [T3] The maximum frequency following the first 
backscattering is given by LOout/^m = (1 + |/3|)/(1 ~ |/3|)^ and following the second by [(1 + |/3|)/(1 - |/3|)]^. One also sees the 
effect of downscattering to a frequency [(1 - |/3|)/(1 + |/3|)]tJin, since some of the charges are flowing outward (in at least one 
hemisphere). The low-frequency spectral slope is generally harder than the Rayleigh-Jeans value of 2, which means that the 
output spectrum has a nearly Rayleigh-Jeans shape at low frequencies when the input spectrum is a single-temperature black 
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Fig. 1 3 . — Left: Output spectrum in a twisted magnetosphere (A0n-s = 1) obtained by integrating eq. 147 1 over magnetic polar angle. Seed photons are injected 
in either the E-mode (solid lines) or the O-mode (dashed lines). Vertical dotted lines mark out the maximum frequency shift in each order of scattering; they are 
separated by factors of (1 + |/3|)/(1-|/3|) in frequency. Right: Angular distribution of the unscattered photons (lOom = i^in) with respect to the magnetic polar 
axis. In both cases the top plots assume a mono-energetic particle distribution (eq. 1361 ) with mean drift speed /3 = 0.5, total width A/3 = 0.1, and unidirectional 
particle flow. The bottom plots assume the same distribution function, but a bidirectional particle flow. 

body. The observation of a flatter spectrum at lower frequencies in some AXP sources jDurant & van Kerkwiikll2006al) therefore 
points to the presence of temperature gradients across the neutron star surface. 

The distribution of number of scatterings n is shown in the left panel of Fig. [T4]for the case of undirectional particle flow (and 
the same simulation parameters as in Fig. [T3] l. Matching the power-law slope with eq. lAli one infers a small average optical 
depth per scattering, r ~ 0.02, for A0n-s = 1 and (3 = 0.5. Both seed photon polarizations produce nearly the same slope. The 
simulation includes enough photons (A^ = 10^) to obtain good statistics for the first ~ 20-25 scatterings. 

The output energy spectrum resulting from the input of a single-temperature black body is plotted in the right panel of Fig. [14] 
Here two methods are used, which give nearly identical results. The continuous line shows the result obtained by drawing the 
seed photon frequencies directly from a blackbody number distribution (eq. P6l ); the dashed line shows the result obtained by 
using a single input frequency (wpeak, eq.|49]l, and then convolving the output spectrum with a blackbody. The main advantage of 
the second method is that the output spectrum is somewhat smoother, and we adopt it throughout this paper. 

The energy spectra obtained with the current code were cross-checked with the output from a Monte Carlo code written 
previously by one of the authors (CT), which assumed a (split)-monopole magnetic field geometry, a mono-energetic particle 
spectrum, and calculated the position of successive scattering surfaces by a different (geometric) method. This previous code 
treated the transport of a photon within a resonant surface in a one-dimensional approximation. The two codes gave identical 
results for the high-energy power-law index, given the same field geometry, angular distribution of optical depth (sin^ 0; eq. 
ITTI ). and particle drift speed. The output spectrum differed only in some details in the shape of the peaks at the first and 
second scattering orders, which were ascribable to the neglect of the spatial curvature of the resonant surface in the 1-D code. In 
particular, a photon that is scattered in a direction tangent to the resonant surface sees an unrealistically high optical depth when 
the surface curvature is neglected. This resulted in the presence of a spurious dip at a frequency (Wout/win) = 1/(1 - in the 
1-D output spectrum. 

Line-of-Sight Dependent Quantities 

In this paper we calculate pulse profiles in different frequency bands (eq. 11611 ). and along specific lines of sight. Shown in the 
left panel of Fig.[T5]is the angular distribution of emission, eq. (ISOl l, for the same simulation shown in Figure[T3] assuming E-mode 
seed photons. These angular distributions are asymmetric with respect to ji^ = 0, which allows for large pulsed fractions. The 
higher photon flux in the south magnetic hemisphere (/i^ < 0) reflects the southward drift of the charges near the magnetic equator 
(where the optical depth is largest). The angular distribution in the two lowest frequency bands also reflects the contribution of 
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Fig. 14. — Left: Distribution of number of scatterings for the simulation shown in the upper panels of Fig. [13] As expected in a multiple scattering process 
with optical depth independent of the order of scattering (eq. IA2I ). the distribution is exponential. Both seed polarizations generate the same slope, although the 
number of unscattered photons is bigger for the O-mode due to its smaller mean optical depth. Right: Energy spectrum for the same simulation parameters and 
E-mode seed photons, but with a di ffere nt handling of the input frequency distribution. The solid line shows the result obtained by drawing frequencies directly 
from a blackbody distribution (eq. 1461 ). whereas the dashed is obtained by convolving the response function of a monochromatic seed spectrum (frequency 
2.?>2k^T\ib/'h; left panel of Fig. [13) with a blackbody distribution. Energy spectra are normalized so that the area of the Planck function is unity. 




Fig. 15. — Left: Angular distribution of outgoing photons (eq. 1501 ) in the five energy bands defined in eq. )6U . Bands III and IV contain most of the photons 
in the non-thermal tail of the spectrum, which is generated by multiple resonant scattering. Right: Pulsed fraction as function of beam angle Ohs^m obtained from 
the angular distributions shown in the left panel. Convergence to 1% in all bands is achieved for Sbeam = ■'r/20 = 9°. 

unscattered seed photons, which peaks along the magnetic poles (see the right panel of Fig.[T3]l. A much larger fraction of the 
photons in bands III and IV have undergone at least one scattering. 

Because only a finite number of photons are evolved in each simulation, a pulse profile can be obtained only by averaging the 
photon flux over some angular cone of half-width ^beam (as in eq. 1521 ). In reality, (?beam — Rks/D ~ 10"'^ for a source of size 
^res ^ 100 km (eq. ||6l) at a distance D ^ 1> kpc. In practice, ^beam is taken to be the largest value for which the pulse profile 
has essentially converged to a constant shape, and faithfully reproduces the result for ^beam — 0. As a measure of convergence, 
we calculated the pulsed fraction as a function of beam width for a few respresentative profiles, and found convergence to within 
~ 1 percent at ^beam = 7r/20 = 9 deg (see the right panel of Fig. [Tsl l. In this calculation the angular grid is taken to have a size 
dji' ~ d(j)' ~ (1 - /Ltbeam)/2 — (^beain/2)^, which yields a large number of angular cells within the beam. 
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